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SUMMARY 

A FORTRAN IV subprogram, WASP, was developed to calculate the thermodynamic 
and transport properties of water and steam. The temperature range is from the triple 
point to 1750 K, and the pressure range is from 0. 1 to 100 MN/m 2 (1 to 1000 bars) for 
the thermodynamic properties and to 50 MN/m 2 (500 bars) for thermal conductivity and 
to 80 MN/m 2 (800 bars) for viscosity. WASP accepts any two of pressure, temperature, 
and density as input conditions. In addition, pressure and either entropy or enthalpy are 
also allowable input variables. This flexibility is especially useful in cycle analysis. 

The properties available in any combination as output include temperature, density, 
pressure, entropy, enthalpy, specific heats (C and C y ), sonic velocity, (3P/3 p) T , 
(3P/9T)p, viscosity, thermal conductivity, surface tension, and the Laplace constant. 

The subroutine structure is modular so that the user can choose only those subrou- 
tines necessary to his calculations. Metastable calculations can also be made bv using 
WASP. 


INTRODUCTION 

Water is inert, inexpensive, and available. It is used for cooling equipment, for 
heating or cooling other fluids, as a modeling fluid, and in many cases as the primary 
test fluid in heat-transfer and fluid dynamics research. 

Printed tables of water and steam properties have been available to the engineer for 
many years, the latest accepted editions being references 1 and 2. Numerous computer 
codes to interpolate these tables using a variety of curve-fit and interpolation techniques 
are available. Many are cumbersome or lack the ability to calculate a consistent set of 
properties for a given point in the fluid surface. Some are designed for specific uses 
and do not include all the properties. A comprehensive, flexible, and internally con- 
sistent computer code for water properties was needed at the Lewis Research Center. 
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In determining the coefficients of equation (1), the temperature data were all ex- 
pressed on the thermodynamic Celsius temperature scale. Since experimental observa- 
tions for water, however, are usually reported in terms of the International Practical 
Scale (the I. P. Scale), a lengthy discussion and a graph of the relation between the 
I. P. Scale and the thermodynamic temperature scale are presented in reference 3. ^ 
The critical constants of reference 3 differ from those presented in references 1 
and 2 as follows: 



Refer- 

References 


ence 3 

1 and 2 

Critical pressure, 

22. 089 

22. 120 

P c , MN/m 2 



Critical tempera- 

374.02 

374. 15 

ture, T , °C 



Critical density, 

0. 317 

0. 31546 

P c , g/'cm 3 




The temperatures in this table are on the I. P. Scale. The critical temperature T of 
reference 3 on the thermodynamic scale is 374. 136° C. C 

The WASP subprogram was developed to be used in fluid-flow and heat-transfer 
calculations. There are independent calls for obtaining any one of the three state vari- 
ables (pressure, density, and temperature) as a function of the other two (see table I, 
OPERATIONS SHEET FOR SUBROUTINE WASP). In addition, temperature and all the 
other properties can be obtained as a function of pressure and enthalpy (or pressure and 
entropy). This option is of considerable value in forced-convection studies and cycle 
analysis. 

While enthalpy, entropy, and specific heats (C and C y ) are available in refer- 
ence 3, the sonic velocity, viscosity, thermal conductivity, and surface tension were 
not computed in reference 3. The sonic velocity, equation (B30), is defined in terms of 
equation (1). The transport properties are discussed in the following section. 

TRANSPORT PROPERTY CALCULATIONS 


The thermal conductivity, viscosity, and surface tension are available in references 

2 

Differences between the current conversion from Celsius to I. P. Scale and that of ref. 3 
are small except at elevated temperatures. These deviations did not warrant our reexamination 
of equation (1). 


4 




1 . Report No. 2. Government Accession No. 

NASA TN D-7391 


4. Title and Subtitle 

WASP - A FLEXIBLE FORTRAN IV COMPUTER CODE FOR 
CALCULATING WATER AND STEAM PROPERTIES 


7. Author(s) 

Robert C. Hendricks, IldikoC. Peller, and Anne K. Baron 


9. Performing Organization Name and Address 

Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio 44135 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D. C. 20546 

15. Supplementary Notes 


3. Recipient's Catalog No. 


5. Report Date 

November 1973 


6. Performing Organization Code 


8. Performing Organization Report No. 

E-7339 


10. Work Unit No. 

502-04 


11. Contract or Grant No. 


13. Type of Report and Period Covered 

Technical Note 

14. Sponsoring Agency Code 


16. Abstract 

A FORTRAN IV subprogram, WASP, was developed to calculate the thermodynamic and trans- 
port properties of water and steam. The temperature range is from the triple point to 1750 K, 
and the pressure range is from 0. 1 to 100 MN/m 2 (1 to 1000 bars) for the thermodynamic prop- 
erties and to 50 MN/m 2 (500 bars) for thermal conductivity and to 80 MN/m 2 (800 bars) for 
viscosity. WASP accepts any two of pressure, temperature, and density as input conditions. 

In addition, pressure and either entropy or enthalpy are also allowable input variables. This 
flexibility is especially useful in cycle analysis. The properties available in any combination 
as output include temperature, density, pressure, entropy, enthalpy, specific heats (C and 
C v ), sonic velocity, (9P/3 p)^,, (9P/3T)^, viscosity, thermal conductivity, surface tension, 
and the Laplace constant. The subroutine structure is modular so that the user can choose 
only those subroutines necessary to his calculations. Metastable calculations can also be made 
by using WASP. 


17. Key Words (Suggested by Author (s)) 

Thermodynamics 

Computers 

Properties of water and steam 


18. Distribution Statement 

Unclassified - unlimited 


19. Security Classif. (of this report) 

Unclassified 


20. Security Classif. (of this page) 

Unclassified 


21. No. of Pages 
118 


* For sale by the National Technical Information Service, Springfield, Virginia 22151 









WASP is a FORTRAN IV subprogram developed for engineering calculations. The 
thermodynamic properties are calculated by using the Helmholtz free-energy equation 
developed by Keyes, Keenan, Hill, and Moore (ref. 3). The transport properties are 
calculated by using curve fits given in references 1 and 2 in regions where these equa- 
tions exist. The authors developed their own approximations based on the tabulated 
values of references 1 and 2 where curve fits were not available. 

The main section of this report is directed to the research-oriented user of WASP. 

It includes a brief discussion of the equations used in calculating thermodynamic and 
transport properties. Comparisons to the International Skeleton Tables and the validity 
of transport calculations are also discussed. A detailed presentation of user instruc- 
tions is included together with a tabular summary for later reference. Detailed informa- 
tion about the computer program and the equations used are included as appendixes. The 
symbols are defined in appendix A; the property equations of WASP are given in appen- 
dix B; the important subroutines of WASP are described in appendix C; the modular de- 
sign of WASP is presented in appendix D; the program listing and flow chart are pre- 
sented in appendix E, the test program output in appendix F, the metastable subroutine 
PMETAS in appendix G, and the thermodynamic relations and derivatives in appendix H. 


THERMODYNAMIC CALCULATIONS 

Keyes, Keenan, Hill, and Moore (ref. 3) fit the available experimental water and 
steam data from the triple point to a pressure of 100 MN/'m 2 and to a temperature of 
about 1750 K, using the fundamental equation 

^ = t|/(T,p) 

= i^ 0 (T) + RT [In p + pQ(p, t)] ( 1) 

where ip is the specific Helmholtz free energy and j = 1000/T. The specific forms of 
^ 0 (T), Q(p,t), and the derivatives of Q(p, t) are presented in appendix B. 

Most investigators (e. g. , refs. 4 to 7), in order to represent their measured values 
as closely as possible, have selected a modified virial equation of state 

P = P(T,p) 


i=l 3=1 


B.(T)p 2i+1 e" cp 


(2) 


2 



where P is the pressure and the coefficients A^T) and Bj(T) are usually polynomials 
in T and T . 

While the derivation of pressure from equation (1) is quite simple, 


P = P 



T 


= p rt(i + pQ + p 2 ^ 


(3) 


the ensuing expanded descriptions for both equations (2) and (3) are quite involved, see 
appendix B and reference 8. 

One should note the influence of the modified forms of the Benedict-Webb-Rubin 
(BWR) equation of state (see refs. 4 to 6) and the more recent work of Bender (ref. 7) 
on equation (1). (Compare the form of Q in eq. (B4) of appendix B with eq. (2).) The 
authors of references 5 and 6 added new exponential terms to the BWR equation of refer- 
ence 4 to account for high-density effects. The technique has been successfully applied 
to several cryogens. More recently, Bender (ref. 7), in addition to these modifications, 
imposed another constraint, namely that the Maxwell Phase Rule must be satisfied; the 
constraint requires that the Gibbs free energy for the saturation liquid and vapor be 
equal. This latter constraint, although not stated explicitly, is implicitly satisfied by 
equation (1) (taken from ref. 3) because the Gibbs free energy of the saturated liquid and 
vapor are ’’virtually identical. ” 

Both equations (1) and (2) have been fit by using a weighted least-squares technique 
which minimizes the residuals in pressure subject to various constraints such as 



T = T 

c 


at the thermodynamic critical point. Reference 3 cites 14 such constraints; usually, the 
number is about one -half as many. However, the advantage of the reference 3 approach 
is that \p as a function of p and T is a fundamental equation and all thermodynamic 
properties are obtained directly from \p and its derivatives. In equation (2), P as a 
function of p and T is a state equation. In determining properties such as enthalpy, 
entropy, and specific heats, the state equation must be differentiated and integrated and 
the associated constants of integration must be determined from other data. * 

*The mathematical form of the derived and integrated equations must be such that they do not 
possess singularities except at the critical point. 
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In determining the coefficients of equation (1), the temperature data were all ex- 
pressed on the thermodynamic Celsius temperature scale. Since experimental observa- 
tions for water, however, are usually reported in terms of the International Practical 
Scale (the I. P. Scale), a lengthy discussion and a graph of the relation between the 
I. P. Scale and the thermodynamic temperature scale are presented in reference 3. ^ 
The critical constants of reference 3 differ from those presented in references 1 
and 2 as follows: 



Refer- 

References 


ence 3 

1 and 2 

Critical pressure, 

22.089 

22. 120 

P c , MN/m 2 



Critical tempera- 

374.02 

374. 15 

ture, T c , °C 



Critical density, 

0. 317 

0. 31546 

P c , g/cm 3 




The temperatures in this table are on the I. P. Scale. The critical temperature T c of 
reference 3 on the thermodynamic scale is 374. 136° C. 

The WASP subprogram was developed to be used in fluid -flow and heat-transfer 
calculations. There are independent calls for obtaining any one of the three state vari- 
ables (pressure, density, and temperature) as a function of the other two (see table I, 
OPERATIONS SHEET FOR SUBROUTINE WASP). In addition, temperature and all the 
other properties can be obtained as a function of pressure and enthalpy (or pressure and 
entropy). This option is of considerable value in forced-convection studies and cycle 
analysis. 

While enthalpy, entropy, and specific heats (C p and C y ) are available in refer- 
ence 3, the sonic velocity, viscosity, thermal conductivity, and surface tension were 
not computed in reference 3. The sonic velocity, equation (B30), is defined in terms of 
equation (1). The transport properties are discussed in the following section. 

TRANSPORT PROPERTY CALCULATIONS 

The thermal conductivity, viscosity, and surface tension are available in references 

2 

Differences between the current conversion from Celsius to I. P. Scale and that of ref. 3 
are small except at elevated temperatures. These deviations did not warrant our reexamination 
of equation (1). 
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1 and 2, with the exception of a small region near the thermodynamic critical point for 
thermal conductivity and a large region (573 to 647. 3 K) for viscosity. For these re- 
gions, the tabulated values of references 1 and 2 represent the output of some interpola- 
tion scheme. However, no technique for predicting transport properties in these regions 
is given in these references. 


Viscosity and Thermal Conductivity 

Reference 8 uses the simple empirical technique of reference 9 to express the vis- 
cosity 77 and thermal conductivity X of several fluids. An excess function of density p 
is added to the dilute gas function of temperature T such that 


n = m + Ar] (4) 

Atj = 7} - = F(p) (5) 

1 U = F(T) at P = 0. 1 MN/m 2 (6) 

and 

X = X x + AX (7) 

AX = X - X x = F(p) (8) 

X 1 = F(T) at P = 0. 1 MN/m 2 (9) 


We used this technique to obtain viscosity and thermal conductivity in the regions where 
no equation existed in references 1 and 2. 

The Atj’s shown in figure 1 were obtained from the tabulated and computed data in 
references 1 and 2. These data were then extrapolated through the region where no fit 
is available (region 4) to predict the viscosity in that region. 

Figure 1 gives a good representation for viscosity over a considerably larger range 
than 573 to 647. 3 K, with deviations from the tabulated data up to 7 percent in some re- 
gions and perhaps 10 percent in the critical region. References 1 and 2 list uncertainties 
at ±4 percent in those regions where the values were not interpolated, and no uncertain- 
ties are given for the interpolated region. The analytical representations of viscosity 
are given in appendix B. 

The AX's used to predict thermal conductivity were also obtained from the tabu- 
lated and computed data of references 1 and 2. In this case the scatter is more acute 
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over a wide range of density because these data are pressure sensitive and are not en- 
tirely represented by excess thermal conductivity as a function of density, see figure 2. 
While the curve of figure 2 represents a wider P, T range than the region where no 
curve fit is available (region V) in figure 2, the AA curve fit is used only in region V. 

Generally, deviations to 8 percent in A -calculations can be found with respect to 
the tabulated data. References 1 and 2 give the deviations as ±4 percent in regions 
where curve fits exist, and no uncertainties are listed for the interpolated region. The 
analytic forms for thermal conductivity are given in appendix B. In region IV (fig. 2) 
the implicit equation for thermal conductivity from reference 2 is used. In region III an 
explicit expression for thermal conductivity from reference 2 is used. These forms 
were adopted over those of reference 1 because of their analytic nature. 

The 0. 1-MN/m 2 thermal conductivity and viscosity output from 700 to 1700 K was 
checked against the results of Svehla (ref. 10). Svehla’ s viscosity is 5 to 10 percent 
higher and his thermal conductivity is 10 to 15 percent higher than those predicted here- 
in. Since the publishing of reference 10, Svehla has found that inclusion of a rotational 
relaxation effect would lower his calculated viscosity perhaps 5 percent, and lower his 
calculated thermal conductivity perhaps 10 to 15 percent. 


Near-Critical Thermal Conductivity 

The anomalous behavior of thermal conductivity in the near-critical region was 
measured by Le Neindre (ref. 11). Sengers (ref. 12) advanced a technique to predict 
the behavior of near -critical thermal conductivity data for carbon dioxide. In refer- 
ence 13, Sengers' technique was modified, extended to several fluids, and compared with 
other methods. The technique used herein is the same as Sengers' technique as pre- 
sented in reference 13, except that the proportionality constant, 3. 05x10 as given by 
Sengers in reference 12, has been increased to 11.6X10’ 5 . 


THERMODYNAMIC AND TRANSPORT PROPERTY PLOTS 

Sample plots of the properties calculated by WASP are found in figures 3 to 14. 2 
The triple temperature scales in these figures are to facilitate checking. Figure 3 rep- 
resents density as a function of temperature for selected isobars, including the critical 
isobar. No irregularities were found. Figure 4(a) represents the pressure as a func- 
tion of temperature for selected isochores. These isochores exhibit distinct curvature 
not only near th e saturation boundaries but in the extended regions as well. The slopes 

3 

Isobars which cross the saturation boundaries are parallel to isotherms. The plotting routine 
simply connects increments in temperature. 
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of the isochores as a function of temperature are shown in figure 4(b), which reveals the 
nonlinear character of most of the isochores of figure 4(a). Figure 5 represents the 
pressure-volume (P-V) plane for selected isotherms. Local P-V regions could be 
mapped by using WASP for preliminary cycle analysis. Figure 6 represents enthalpy as 
a function of temperature for selected isobars, including the critical isobar. Figure 7 
represents the temperature -entropy (T-S) diagram for selected isobars. Again, local 
T-S regions could be mapped by using WASP for preliminary cycle analysis. 

Figure 8 represents specific heat at constant pressure and figure 9(a) represents 
specific heat at constant volume for selected isobars, including the critical isobar. Note 
the peaking effects in C y along the critical isobar; this indicates a discontinuity in C , 
as well as in C , along the critical isobar. This behavior agrees with the most recent 
thinking on C„ in the critical region; namely, that C, r possesses a weak logarithmic 

W V V 

infinity (i.e., C a |T ■ T |" ' , as discussed in ref. 14). In figure 9(b), the isen- 

tropic expansion coefficient (y = Cp/C y ) is given as a function of temperature for se- 
lected isobars, including the critical isobar. Note that while Cp and C y tend toward 
a discontinuity, so also does the ratio. The reason is that C possesses a very strong 

( \ -h p \ 

"infinity” at the critical point ( C ^ T - T c , where 1. 2 < j8j < 1. 35 j and C y a 

rather weak "infinity. " Consequently, in the critical region, y diverges approximately 

- 0 , 

as T - T c where 1. 1 < 0 2 < 1. 3. 

Figure 10 gives sonic velocity as a function of temperature for selected isobars. 
Sonic velocity tends toward zero, or at least a minimum, in the critical region C , 
which also concurs with recent thinking because 


c c - 


T /9P 


Since C y diverges in a weak manner and (9P/9T)p is nonzero and finite, C c will ap- 
proach zero in a weak manner. 

Figure 11 represents a plot of viscosity as a function of temperature along selected 
isobars, including the critical isobar. The discontinuities of this surface at 573 and 
647 K are caused by the empirical nature of the curve fit (fig. 1), as discussed previ- 
ously. This discontinuity is on the average less than 4 percent, which is the same order 
of accuracy as presented in references 1 and 2. 

Figure 12(a) represents thermal conductivity as a function of temperature along se- 
lected isobars. In this report, an attempt has been made to include the anomalous be- 
havior of thermal conductivity in the near-critical regions based on references 12, 13, 
15, and 16. The behavior of the near-critical thermal conductivity is shown in fig- 
ure 12(b). In order for the user to obtain these values, he must add EXCESSK, which 
represents the anomalous part of the thermal conductivity, to the normally computed 
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value of thermal conductivity. See USER'S GUIDE TO WASP, the statement COMMON, 
PROPTY/ .... Generally, the plots of ?j and A exhibit some irregularities where 
the various predicting techniques overlap; however, for most applications the values re- 
turned are within acceptable tolerance. 

Figure 13 is a plot of surface tension a and Laplace constant as a function of tem- 
perature. Metastable conditions near both the liquid and vapor are often required in a 
system analysis and can be calculated by equation (1). A special subroutine is included 
in appendix G which when used with WASP will give metastable properties; sample plots 
are shown as figure 14. 


Other Thermodynamic Functions 

WASP provides sufficient PVT and derived property data for most users; however, 
if other functions are required, the user may calculate these by using the partial deriva- 
tives (3P/3p) T and (0P/3T)^, along with the other results from WASP. Appendix H is 
provided to give the user a handy reference to the so-called Bridgeman Tables which 
list most of the interrelations between thermodynamic variables. (See pp. 36 and 64. ) 


Comparison Plots 

Of utmost importance is how well equation (1), as used in WASP, agrees with the 
International Skeleton Tables for steam and water (refs. 1 and 2). Figures 15 to 17 were 
obtained by running all the data points listed in table 1. 2 of reference 1 (see also table 4b 
of ref. 2) on each of the three input options (T, P), (p, T), and (p, P) in WASP. Each of 
the figures is discussed, but the careful reader should note that discrepancies exist in 
the specific volumes presented in these two references. The authors noted four obvious 
errors in reference 2 by comparing references 1 (table 1. 2) and 2 (table 4b). Refer- 
ence 1 was assumed to be correct. Other discrepancies occur in reference 2, for in- 
stance in the specific heat C . No attempt was made to track all these errors, and the 
reader should use reference 2 cautiously. 

Figure 15 represents the percent relative error in density, [(p^ie - P C alc^ 
p table-l X 100 ’ as a function of density. With the exception of three points, all the 
values are within +0. 25 percent and -0. 50 percent, and generally have an error of less 
than 1 part in 3000. The solution for density is iterative, and perhaps the error could be 
reduced somewhat by a tightening of the convergence criteria. This is not recommended 
for two reasons: (1) It will require a great increase in computer time, and (2) errors in 
these printed tables have been noted. The tolerance should be quite satisfactory to all 
but the most critical user. 
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Figure 16 represents the percent relative error in pressure, [.P ta jj le - p ca ^ c )/ 
^table-1 X as a ^ unc ^ on °f pressure. In all cases the calculated pressures are 
within +3. 0 percent and -2. 0 percent of the tabulated value; most points lie within +0. 25 
percent of the tabulated value. The prediction of pressure at high density (low tempera- 
tures) using a fundamental equation or a state equation is quite difficult. These pres- 
sure errors are all within accepted tolerances. 

Figure 17 represents the percent relative error in temperature, Citable “ T calc^ 
^table-1 X as a f unc ti° n of temperature. With the exception of about a dozen points, 
the predicted temperatures are within +0. 25 percent and -0. 4 percent and generally lie 
within ±0. 1 percent. 

Usually, temperature and density are predictable because of the manner in which 
the data were acquired; however, pressure is always difficult to calculate. With these 
basic guidelines in mind and figures 15 to 17, it can be said that the equation gives a 
faithful representation of the International Skeleton Tables (refs. 1 and 2). 


USER'S GUIDE TO WASP 

The user with limited programming experience should have no difficulty in following 
the operating instructions for WASP. After gaining a little experience with WASP, the 
only references needed are table I (the operations sheet) and table II (the units specifica- 
tion). 


How WASP Handles Input/Output 

WASP is a group of subroutines designed to be used as a subprogram with the user's 
program. Standard communication between the user's program and WASP is achieved 
by the following two FORTRAN statements, which contain the symbols representing the 
input/output parameters and options: 

COMMON/PROPTY/KU, DL, DV, HL, HV, etc. 

CALL WASP (KS, KP, T, P, D, H, KR) 

See table I and appendix A. 

Three requirements must be fulfilled for a successful call to WASP: 

(1) The cards for COMMON/PROPTY/KU, DL, DV, etc. , must be included in the 
user’s main program or the subroutine that calls WASP. The WASP subprogram deck 
must be correctly loaded with the user’s program as shown in table in. The variables 
MU, MUL, MUV, K, KL, and KV must be declared REAL. (K cannot be used as an 
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index for a subscripted variable. ) However, the user can change the names of these 
variables in COMMON/PROPTY/KU, DL, DV, etc. , . . . . 

(2) The units system for input/output must be correctly specified. KU is an input 
control specified in the COMMON/PROPTY/KU, DL, DV, etc. , which must be set such 
that 1 ^ KU < 5. KU identifies the units system for input/output, and KU is never al- 
tered by WASP. Therefore, unless the user switches from one system to another, he 
need set his value for KU only once, before any calls to subroutine WASP. 

There are three specified units options described in table II. The option KU=I is 
the internal program units system. The other two options are commonly used in en- 
gineering calculations. If the user does not wish to use one of these options, he can 
specify any desired units system for KU=4 and KU=5, provided the conversion factors 
for this system are stored by the user as directed in table n. 

(3) The controls KR, KS, and KP, which tell WASP what variables are to be used 
as input and what properties are requested for output, must be correctly initialized in 
the call statement for subroutine WASP. The corresponding input variables in the call 
statement and COMMON/PROPTY/ . . . must also be correctly initialized. 

KS and KR are controls that determine which of the variables T, P, D, H, or S or 
combinations thereof are needed as thermodynamic input. KP is an input control which 
specifies which properties are sought as output. 

KR is also an output variable since it gives the correct region number for the vari- 
ables in a specific call to WASP, as shown in the sketch in table I. Depending on the 
input for KS and KP, the other possible output variables are T, P, D, H, and all of 
COMMON/PROPTY/ except the control KU. 

As mentioned above, KR is both input and output and must be reset before each call 
to WASP. The input options are 

(a) KR=0 when user wishes WASP to determine a value for KR 

(b) KR=1 when user wishes saturation conditions 4 
The output for KR will be 

(a) KR=1 for saturation 

(b) KR=2 for liquid 

(c) KR=3 for vapor 

KS specifies which variables are to be used as input for a call to subroutine WASP. 
(In the remaining discussion on WASP input/output, the input variables are assumed to 
be in user's units specified by KU. Output is always returned in the KU system of units. ) 

^Saturation or coexistence conditions exist on the PVT surface when pressure is a function of 
only temperature and the liquid and vapor states both exist at that pressure. Thus when KR=1, 
two outputs for each property are available in COMMON/PROPTY/ and only one independent 
variable is required for some input options, as shown in the KS-KR input/output chart. 
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The following table shows the input and output for all KS, KR combinations: 


Thermodynamic 
region specifica- 
tion, 

KR 

State relation specification, KS 

1 

2 

3 

4 

5 

Input 

0 

T and P 

T and D 

P and D 

P and H 

P and S 

1 

T or P a 

T 

P 

P 

P 

Output 

1 

T or P a , 

P 

T 

T, DL, DV 

T, DL, DV 


DL and DV 





2 

D 

P 

T 

D and T 

D and T 

3 

D 

P 

T 

D and T 

D and T 


If T is the desired input, set P = 0. 0 prior to the call and vice versa. Then 
WASP will return the correct saturation value for the 0. 0 input. If both T 
and P have an input value t 0. 0, WASP uses T but will not alter P input. 


an input control that specifics which derived and transport properties are re 
quested by the user. It is the sum of the individual KP options and is described in ta- 
ble I. This binary sum allows WASP to uniquely identify any combination of requests. 
The following table shows the output locations for the specific KR and KP combinations: 


Value added 
to KP input 

Output for 
KR=2 or 3 

Output for KR=1 

Name of calculated property 

Liquid 

Vapor 


0 


--- 

--- 

— 

None requested 

1 


H 

HL 

HV 

Enthalpy 

2 


S 

SL 

sv 

Entropy 


r 

CP 

CPL 

CPV 

Specific heat at constant pressure 

4 i 


cv 

CVL 

CVV 

Specific heat at constant volume 



GAMMA 

GAMMAL 

GAMMAV 

Specific-heat ratio 


L 

C 

CL 

CVP 

Sonic velocity 

8 


MU 

MUL 

MUV 

Viscosity 

16 


K 

KL 

KV 

Thermal conductivity 

1 

r 

SIGMA 

SIGMA 

--- 

Surface tension of the liquid as a 

32 < 





function of temperature 

1 


ALC 

ALC 

— 

Laplace constant as a function of 


1 




temperature 
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Troubleshooting for User Errors 


After experience with WASP, we have found that several common errors are easily 
detected and corrected. 

(1) Failure to set 1 < KU < 5 will cause a "division by 0. 0” and/or no valid 
answers. Set KU to its proper value. 

(2) Failure to set 1 < KS < 5 will most likely cause a halt to the program because 
of an execution error. The branching on KS in subroutine WASP is a computed "GO 
TO. " Simply set KS to its proper value. 

(3) Failure to set KP will return enthalpy if KP is odd and no derived properties if 
KP is even. 

(4) If a wrong value is entered for KR, it is treated as KR=0. If a user enters 
KR=1 when he does not want saturation properties, he will get them anyway for T < T c 
and otherwise will get a wrong answer. 

(5) If any T, P, D, H, or S is entered incorrectly, that value will be used and the 
answer will be wrong. 

(6) If the COMMON/PROPTY/ is duplicated incorrectly, there are a variety of pos- 
sible errors, almost all serious. 

Other small problems may be encountered if WASP is modified for different com- 
pilers or computers. The FORTRAN IV coding in WASP is machine independent except 
for a few Hollerith format statements which can be easily changed. The reader who 
needs more detailed information should read the appendixes. 


Additional Information 

The approximate core storage for the complete WASP program is (14650) 8 = (6568) 1Q 
locations. 

The time estimates were obtained by running an average of 100 calls over the entire 
PVT range for each option indicated. The shortest call was for pressure, KS=2, at an 
average of 4 milliseconds per call for T > T and 17 milliseconds for T < T . The 

^ C 

call for density, KS = 1, varied from 17 to 40 milliseconds for all regions, with the 
greatest time being consumed in the near -critical region. The call for temperature, 
KS=3, varied from 11 to 70 milliseconds per call, with the least time used when P > P 

c 

and the most time used in the near-critical region. The call for density and all the de- 
rived properties, KS=1 and KP=63, varied from 38 to 120 milliseconds per call depend- 
ing on the density call and the regions for the transport properties. 

The P, H and P, S calls, KS=4 and KS=5, each required from 300 to 800 milliseconds 
per call, with the greatest time in the near -critical region. These results are sum- 
marized as follows: 
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State relation specification, KS 

1 

2 

3 

4 

5 

1 

Thermodynamic and transport properties specification, KP 

0 

0 

0 

0 

0 

63 



Time, 

msec/call 



17 to 40 

4 to 17 

11 to 70 

300 to 800 

300 to 800 

38 to 120 


Problems Previously Encountered When Converting to Non-IBM Machines 
or Different FORTRAN IV - FORTRAN V Compilers 

The problems encountered in converting to different equipment are as follows: 

(1) IBM 360 users should run in double precision by inserting implicit REAL+8 
(A-H, O-Z) and REAL* 8 MU, MUL, MUV, K, KL, KV in subroutine WASP and implicit 
REAL*8 (A-H, O-Z) in all other subroutines. Change COMMON/PROPTY/KU, KZ, DL, 
DV, etc. , for proper alinement. 

(2) Data statements are found in subroutines BLOCK DATA, THERM, VISC, and 
SURF. Many compilers differ in formating data statements. 

(3) The multiple-entry routine (CHECK, TCHECK, PCHECK, DC HECK) has an en- 
try point, DCHECK, whose input vector (KU, D) does not correspond in kind and number 
with the other entry points (KU, KR, T) or (KU, KR, P). To our knowledge this has 
caused a problem on only one compiler, a FORTRAN IV for a CDC 3800. It was easily 
remedied by an equivalence statement. 

The authors adapted the code to fit the following compilers and machines: UNIVAC 
1108, CDC 3600, CDC 6600, IBM 360/67TSS, and IBM 7094-7044 DCS. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, April 26, 1973, 

502-04. 
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APPENDIX A 


SYMBOLS 


Mathemat- 

ical 

symbol 

A.. 

ij 


Di 


D *i 

E 

H 


FORTRAN 

symbol'* 


A(I, J) 

ALC 

C 

CL 

CVP 

SICl^ 


SIC5J 


CP 

specific 

CPL 

specific 

CPV 

specific 

CV 

specific 

CVL 

specific 

CVV 

CPSfl 

specific 


CPS7_ 

E 

H 

HL 


Definition 

coefficients of terms in Q-function (see table IV) 

Laplace constant 

sonic velocity, cm/sec 

sonic velocity of saturated liquid, cm/sec 

sonic velocity of saturated gas, cm/sec 

coefficients of terms in function (see table IV) 

specific heat at constant pressure, J/(g)(K) 

at C of saturated liquid, J/(g)(K) 

r 

at C of saturated vapor, J/(g)(K) 

r 


V 

'v 


coefficients for vapor pressure curve (see table IV) 


E = 4. 8 
enthalpy, J/g 

enthalpy of saturated liquid, J /g 


^Symbols used in each individual subroutine are identified in that subroutine (see appendix C). 
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HV 

enthalpy of saturated vapor, J/g 


KP 

thermodynamic and transport properties specification 


KR 

thermodynamic region specification 


KS 

state relation specification 


KU 

units specification 

p 

P 

2 

pressure, MN/m 

Q 

QCALC 

data-fitting function 

R 

R 

specific gas constant for water, 0. 46151 J/(g)(K) 

S 

S 

entropy, J/(g)(K) 


SL 

entropy of saturated liquid, J/(g)(K) 


SV 

entropy of saturated vapor, J/(g)(K) 

T 

T 

temperature, K 

u 


internal energy, J/g 

y 

GAMMA 

ratio of specific heats, Cp/C y 


GAMMAL 

ratio of specific heats for saturated liquid 


GAMMAV 

ratio of specific heats for saturated vapor 

V 

MU 

dynamic viscosity, g/(cm)(sec) 


MUL 

dynamic viscosity of saturated liquid, g/(cm)(sec) 


MUV 

dynamic viscosity of saturated vapor, g/(cm)(sec) 

A 

K 

thermal conductivity, W/(cm)(K) 


KL 

thermal conductivity of saturated liquid, W/(cm)(K) 


KV 

thermal conductivity of saturated vapor, W/(cm)(K) 

p 

D 

9 

density, g/cm 

Pl 

DL 

9 

density of saturated liquid, g/cm 

p v 

DV 

9 

density of saturated vapor, g/cm 

(>* 

RHOA 

constant used in Q-function, p„ = 0. 634 


RHOB 

constant used in Q-function, = 1. 0 

a 

SIGMA 

surface tension, dyne/cm 

1000 

T = 

TAU 

temperature parameter, 


T 
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T a 

TAUA 

T c 

TAUC 

* 

PSI 

^0 

PSIO 

d^o 


PSIT 

dT 


QTD 

\B t / 

p 

(f) 

QDT 



/ cV 

Q2T2D 

U 2 y 


/d 2 Q\ 

P 

Q2D2T 

V/ 


a 2 Q 

T 

Q2DT 

dr dp 

J 


constant used in Q-function, j ^2.5 

1000 divided by critical temperature expressed in kelvin 

Helmholtz free energy, J/g 

reference function, J/g 


partial derivatives used in evaluating \j/ and its derivatives 


16 



APPENDIX B 


PROPERTY EQUATIONS OF WASP 

The equations used in WASP were those taken from Keyes, Keenan, Hill, and Moore 
(ref. 3), Schmidt (ref. 2) and the ASME Steam Tables (ref. 1) and those developed by the 
authors. 


FUNDAMENTAL EQUATION 

The basic equation of WASP expresses the Helmholtz free energy in terms of p 
and T, 


'P = 'Pip, T) 


(Bl) 


whereas the equation of state is usually expressed as 


P = P(p, T) 


(B2) 


Equation (Bl) is complete inasmuch as the required thermodynamic functions are deriva- 
tives of \f/ and undetermined constants and/or functions are not required. For example, 
specific heat at "zero” pressure C , which is a function of temperature, is not re- 

p o 

quired in the i//-form; however, in the P-form (eq. (B2)), C is required to obtain 

Pq 

entropy, enthalpy, and specific heats. 

When equation (Bl) is expanded, i p becomes 


'P = ^ 0 (T) + RT[ln p + pQ(p, T )] 


(B3) 


where 


8 

Q = Y A il^ p ■ p a^ 1 + e EP ( A 9, 1 + A 10, l p > 

i=l 


+ (t " T c ) \ Y (t ' T a ^~ 2 

j=2 


8 

Y A ij^ p " Pb^ 1 + e EP ( A 9j + A 10j p ^ 

i=l 


(B4) 
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^ 0 (T) =C l + C 2 T + CgT 2 + (C 4 + C 5 T)ln T 

1000 „-l 
7 ~ K 

T 


and the constants C^, 


P a = 0. 634 g/cm 3 
P^ = I- 0 g/cm 3 


T a = 2. 5 K- 1 
T c = 1 . 544912 K' 1 


> 


E = 4. 8 cm 3 /g 


R = 0. 46151 J/(g)(K) 

J 


C & and Ajj are given in table IV. 


(B5) 

(B6) 


(B7) 


DERIVATIVES OF Q 


The derivatives of Q are required to evaluate any of the thermodynamic properties: 



L V p ■ Pi/’ 1 


i=l 


+ e 


-Ep 


(A 9j + 


A 10jP)| 


+ (T - T c ) 




8 


X A ij( p ' Pb) 1 1 + e EP(A 9j 


+ A 



(B8) 



8 

^ (i - 1)A n(p - P a ) + e p [j E(A q ^ + A 1q x p) + A lf) ^ 
i=2 ’ 


+ (t " T c ) 



8 

Z (i - 1)A ij ( p - pb* 1 ' 2 


i=2 


+ e 


-Ep 



+ A 1Qj p)E 
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9 2 Q = 
dp 3 t 


i <' - 'aHt « - « A L1<0 - »b> 

j=2 Li=2 


2 + e Ep (- EA 9j ~ E P A iOj + A lOj^ 


+ (t ‘ t cM Yj G ' 2 Hr - T a ) j ' 3 

U=3 


8 

Z <* - ^(P - P b) 1 ' 2 + e ' EP (- EA 9j * E A 10j P + A 10j ) 
i=2 

(BIO) 


= £ A u (i - l)(i - 2 )(p - p /' 3 + e* Ep {[- EA 0) j + A 10; j(2 - Ep)] (- E)} 


VP/ T i=3 


+ (t - t c )( £ (t - r a ) j - 2 (|; A ij(i - D(i - 2)(P - Pb) 1 ' 3 + e- E P[- EA 9j + A 1Qj (2 - Ep)](- 

J-2 b=3 


E > 


(B11) 


i-Q) =2 ?(j - 2 )( t - r a ) j ' 3 
^r 2 / p U=3 




E A ij (p " p b )1 " 1 + e ’ EP(A 9j + A 10j^ 

i=l 


\J 


+ (T - T c ) 


r 7 

^ (j - 2)(j - 3 )(t - r a )"* 
3=4 


i-4 


E V P " V 1 ” 1 + e " EP(A 9j + A 10j^ 
i=l J 

(B12) 


THERMODYNAMIC PROPERTIES 


The derivatives of i f/ give all the functions necessary to obtain the thermodynamic 
properties. 


Pressure and Its Derivatives 



(B13) 
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(B14) 



Enthalpy and Its Derivatives 

h = u + — 

p 

H = ra(»T)l + p 

L J p P 


H = 

+ 1000 

_ T ^o 

, 1000 R 

[i + pQ + p 2 / i9^ 


' 3r> p 

dT 

T 

L 


where the first term of equation (B19) is the internal energy u 




(B15) 


(B16) 


(B17) 


(B18) 


(B19) 


(B20; 
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Constant volume: 


Constant pressure: 



Specific Heats 



(B25) 


(B26) 



(B27) 
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’’Isentropic” expansion coefficient: 


,3 

(B28) 

C v 


Sonic Velocity 


c 2 = m 
V dp /g 

(B29) 

C 2 = y(— 1 

(B30) 

\ 0P/rp 



Vapor Pressure Curve 

iog 10 p - (i ♦ dj) i5 >’ + tB31 

i=3 

where the original data are in bars and °C whereas pressure and temperature in the pro 
gram are In MN/m 2 and K, hence the forms (1 + Dj) and (T - 273. 15). 


TRANSPORT PROPERTIES 


The transport property equations are not as concisely defined as the fundamental 
equation. The transport maps for viscosity and thermal conductivity are broken in o 
several regions, as shown in sketches a and b, respectively, and individual curve h.s 
are presented for each. Also, several regions are void of description as they exist 

references 1 and 2. 
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^Region 1 



Temperature, T, K 


(a) 

Viscosity 


Atmospheric pressure. - For P = 0. 1 MN/m 2 and 373. 15 K < T < 973. 15 


K, 


^1 " 


b iv b r b3 


xlO 


-6 


R egion 1. - For P gat < P < 80 MN/m 2 and 273. 15 K < T < 573. 15 K, 


(B32) 


rj = 10" 6 a. 1 


* 


x 10 L 


(T/ T c )-a 3 


Regionj. - For 0. 1 MN/m 2 < P < „„ 373 . 15 K < T < 573. 15 K, 


t? = ^xio 6 - io £- 


1 %- c * 


x 10' 


(B33) 


(B34) 


- egl0n 3 ~ ' For °- 1 MN/m 2 < P < 80 MN/m 2 and 648. 15 K < T < 1073. 15 


K, 


V = 


xl ° 


+ d J-&- j + dj£- 


xio 


-6 


(B35) 
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Region 4. - Tabulated values of viscosity in region 4, as well as computed values of 
viscosity at equivalent densities, were plotted as per figure 1. The resulting curve 
gives an accurate representation of these data, with the exception of those values along 
the saturation locus in the near-critical region. As can be seen, deviations of up to 
7 percent, and 10 percent at the critical point, are common. 


k = 1 for p/p £ 4/3 

f 

k = 2 for p/p c > 4/3 J 


(B36) 


n = n i 


10 


0.0192 


(B37) 


where 

Y = ^SkP^ 4 + ^4k x3 + ^k* + C 2k X + C lk 
X V log 10 (£) 


The following coefficients are used in the viscosity equations. 


a x = 241.4 

a 2 = 0. 3828209486 

a 3 = 0. 2162830218 

a 4 = 0. 1498693949 

a c = 0. 4711880117 
o 

b 4 = 263. 4511 
b 2 = 0.4219836243 
b 3 = 80. 4 


Cl = 586. 1198738 
c 2 - 1204.753943 
c 3 - 0. 4219836243 

d x = 111. 3564669 
d 2 = 67. 32080129 
d 3 = 3.205147019 


(B38) 

(B39) 
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C lk = -6. 4556581 


C 2k = 1. 3949436 
C 3k = 0. 30259083 
C 4k = 0. 10960682 
C 5k = 0. 015230031 

For k = 2, 

C l k = -6.4608381 
C 2k = 1.6163321 
C 3k = 0.07097705 
C 4k = -13. 938 
C 5k = 30. 119832 

Thermal Conductivity 

Atmo spheric pressure . - For P - 0. 1 MN/m 2 and 373. 15 K < T < 973. 15 
= (17. 6 + 0. 0587 t + 1. 04xl0” 4 t 2 - 4. 51x10“® t 3 ) xlO" 5 

where 

t = T - 273. 15 

R egion I. - For P gat < P < 50. 0 MN/m 2 and 273. 15 K < T < 623. 15 K, 



(B40) 


(B41) 


(B42) 
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where 


S 


1 " 



s 


2 " 



S 


3 " 



(B43) 


(B44) 


(B45) 


Reglon n . - Far the following ranges of pressure (In MN/m 2 ) and temperature 
(in K): 

0. 1 < P 2= 17. 5 and T gat < T < 973. 15 
17. 5 < P 22. 5 and 673. 15 < T < 973. 15 
22. 5 < P < 27. 5 and 698. 15 < T < 973. 15 

27. 5 < P £ 35. 0 and 723. 15 < T < 973. 15 

35. 0 < P £ 45. 0 and 773. 15 < T <- 973. 15 

45. 0 < P ^ 50. 0 and 823. 15 < T < 973. 15 


thermal conductivity is 


X = 


,14 


-5 2\ 2. 1482X10 2 

X x + (103. 51 + 0. 4198 t - 2. 771x10 t )p + — P 


xlO" 5 (B46) 


where 


t = T - 273. 15 


(B47) 
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and 


Region III. - If the (P, T) is not in region II (see eq. (B54)) but P < 50 MN/m 
373. 15 K < T < 973. 15 K, then the following should be used: 



C = 



c 33 


Region IV . - If the (P, T) is not in region III (see eq. (B54)) but P < 50 MN/m 
T > 623. 15 K, then the following should be used: 



where 


k = 100 X 


The solution for X is iterative. And 



(B48) 

(B49) 


(B50) 


(B51) 

and 


(B52) 


(B53) 
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(B54) 



is the boundary of region III-IV. 

Region V . - In region V, tabulated values of thermal conductivity as well as com- 
puted values of thermal conductivity at equivalent densities were plotted as per figure 2. 
The resulting curve gives a good representation of the tabulated values, except along 
the saturation locus. However, deviations up to 8 percent, and 10 percent near critical, 
can be expected as listed in the table. These tabulated data and this curve fit do not in- 
clude the anomalous behavior of thermal conductivity in the near -critical region. 


k = 1 for < 2. 5 
Pc 


k = 2 for > 2. 5 

P c 


X = X ^ + 10 


(B55) 


(B56) 


(B57) 


where 


C 5k* 


4 + C 4k X 3 + C 3k X + C 2k X 


+ C 


lk 


(B58) 


and 


X = log 



(B59) 
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The constants used in equation (B58) are as follows: 
For k = 1, 

C lk = -0. 5786154 
C 2k = 1- 4574646404 
C 3k = 0. 17006978 
C 4k = 0. 1334805 
C 5k = 0.032783991 

For k = 2, 


C lk = -0. 70859254 
C 2k = 0- 94131399 
C 3k = 0. 064264434 
C 4k = 1. 85363188 
C 5k = 1. 98065901 

The following coefficients are used in the equations for thermal conductivity: 

a Q = -0. 92247 b Q = -0. 20954276 

a l = 6- 728934102 bj = 1. 320227345 

a 2 = -10. 11230521 b 2 - -2. 485904388 


a 3 = 6. 996953832 

b 3 = 1. 517081933 

a 4 = -2. 31606251 

b 31 = 6. 637426916X10 5 

a 31 = 0. 01012472978 

b 32 = 1. 388806409 

a 32 = 0. 05141900883 

b 4Q = 1. 514476538 

a 4Q = 1. 365350409 

b 41 - -19. 58487269 

a 41 = -4.802941449 

b 42 = 113.6782784 



a 42 = 23.60292291 
a 43 = -51. 44066584 
a 44 = 38. 86072609 
a 45 = 33. 47617334 
a 46 = -101. 0369288 
a 47 = 101. 2258396 
a 4g = -45.69066893 

c Q = 0. 08104183147 
c 4 = -0. 4513858027 
c 2 = 0. 8057261332 
c 3 = -0. 4668315566 
c gl = 3. 388557894X10 5 

c 32 = ® 

c 33 = 0. 206 

c 4Q = 1.017179024 


b 43 = -327.0035653 
b 44 = 397. 3645617 
b 45 = 96. 82365169 
b 46 = -703. 0682926 
b 47 = 542.9942625 
b 4g = -85.66878481 

d 31 = 2. 100200454x10 6 

d 32 = 23. 94 

d 33 = 3. 458 

d 34 = 13.6323539 

d 35 = 0.0136 

d 36 = 7.8526X10" 3 

e 4 - 50.60225796 
e 2 = -105.6677634 
e 3 = 55.96905687 


Thermal Conductivity - Anomalous Region 


The Senger technique (ref. 12) as modified in reference 13 and again herein is cal- 
culated for 0. 4 < p/p c — 1-6. Let 
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and \p represent the nonanomalous or frozen thermal conductivity. For 


A - Ap - 


11. 6xl0" 5 



l-.fi- 

Pn 


II. 71 


For X 13 > 3, 


For 0. 4 < < 3, 



where 


x o - 1 

0 

C = iog 10 


and 


a Q = -0. 17384732 
aj = 0. 82350372 
a 2 = -1. 55213983 
a 3 = -0. 12626138 
a 4 = 2. 83922425 


< 0.4, 

(B60) 

(B61) 

(B62) 

(B63) 
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Surface Tension and the Laplace Constant 


Surface tension is given by 


a l( T - T c> 2 
1 - 0. 83(T - T ) 



where 


The Laplace constant is 


a x = 0. 1160936807 
a 2 = 1. 121404688X10' 3 
a 3 = 5. 75280518X10' 6 
a 4 = 1. 28627465X10" 8 
a 5 = 1. 149719240X10" 11 


L = 


{ 


g (P L - P v ) 


where g is the local acceleration. If g is the acceleration of gravity, 

g = 980.665 cm/sec^ 


(B64) 


(B65) 


(B66) 
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APPENDIX C 


description of important subroutines in wasp 


This appendix includes a discussion of the input/output and important features of 

a e7;;:“ eS “ h WASP ' The meth0d 0f SOlUUOn ^ f ° r «“ «*■“« I. indi- 

TMN ,v a r '° equat ‘° ns presented in Wdix B - The FOR- 

RAN IV variables mentioned correspond to the program listing in appendix E. The 

s orter subroutines not included in appendix C are completely described by comments in 

bll I Tn t ‘ X E ' SubrOUtine WASP has been described in the main text, in ta- 
les and II; hence, the reader is assumed to be familiar with subroutine WASP 


MATHEMATICAL ROUTINES 


The mathematical routines are as follows- 

(1) Function SOLVE (XI, F, DF): This routine performs a Newton-Raphson itera- 
tion given the initial estimate XI, the function F, and the derivative function DF The 
convergence is determined when |(X - X )A N | < TOL. The value of TOL is 1 E-5 

or iterations 1 to 40, i. E-4 for 41 60, 1. E-3 L 61 to 80, and 1. E-2 for 81 to 100 

In al l eases studied the convergence was usually obtained In fewer than 40 iterations ' 
For the exceptions, usually in the near-critical region of the PVT surface, the values 
re urned with the increased tolerance are the best obtainable using equation (B3) The 

maximum number of iterations i<? 100 ^ v ine 

number is reached. ’ a PP ro P r ^te message is written if this 

FUNC 2> XD r °Th ineS f°° T (X0, X2 ’ F0FX ’ FUNC ’ X1) a " d ROOTX (X0 > X2 - FOFX - 

essarv for ihe d^h, ! "r" 1 " 5 ‘ denUCal “'"P' t0r name ' The d “Plication is nec- 
essary for the double iterations in the solutions for temperature and density given ores 

sure and enthalpy (KS-4) or pressure and entropy (KS=5) as input. (See also table I ) 

functmn 6 “ IT T 0d ‘ S “ m ° difled haU - i ” lerval “ a -» technique for a monotonic 
on, FUNC, with a root between X0 and X2 such that FUNC(Xl) = FOFX where XI 

is varie n dT„ e t r he e sr ed ' Th6 ' , “ mber ltera “ 0nS d ° eS " 0t eXCeed 100 ' and the France 
varied in the same manner as in function SOLVE. In addition, both the root and the 

function value FUNC(X1, must meet a tolerance. While the tolerance on X, is TOL the 

0 'Tlnn ° n FUNC(X1) ,s Id* TOL. Error messages are written when the iterations 
reach 100 or when there is no solution in the interval X0 to X2. 
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Q -FUNCTIONS 


These routines use D and TAU in program units of KU=1. Entry points with TAU as 
input indicate an iteration where D is known, while entry points with E and TAU as input 
are used in solving for D and in calculating all derived properties. 

The Q-functions are as follows: 

(1) Function QMUST(D) calculates summation terms involving D needed by other 
Q-functions and stores them in COMMON/ QAUX/ and /QSI/. ENTRY QMUST2(TAU) 
calculates summation terms involving TAU and stores them in /QAUX/. 

(2) Function QCALC(TAU) calculates equation (B4). 

(3) Function QTD(TAU) calculates equation (B8). 

(4) Function ODTA(TAU) and ENTRY QDT(D, TAU) calculate equation (B9). 

(5) Function Q2DTA(TAU) and ENTRY Q2DT(D, TAU) calculate equation (BIO). 

(6) Function ©2D2TA(TAU) and ENTRY 02D2T(D, TAU) calculate equation (Bll). 

(7) Function Q2T2D(TAU) calculates equation (B12). 

FUNCTION CHECK 

Function CHECK includes 

(1) ENTRY TCHECK (KU, KR, T) 

(2) ENTRY PCHECK (KU, KR, P) 

(3) ENTRY DC HECK (KU, D) 

These entry points convert the variables from the user's units to the program's 
units, represented by KU=1, and check for out-of-range variables. Appropriate mes- 
sages are written for any out-of-range input, but the calculation is allowed to continue. 

The following subroutines use the mathematical routines, the Q-functions, function 
CHECK, and the subroutines listed with each in table V. The use of these subroutines 
is determined by the KS and KP options (see table I) and are called by subroutine WASP. 
If a user wants to use only a few of these subroutines, he can disassemble the WASP 
program by following the instructions in appendix D and the discussion for the routine of 
interest. Subroutine WASP uses the temperature parameter TAU (in user's units) for 
input to the subroutines. All the derived thermodynamic property and transport property 
subroutines assume that TAU, P, and D have been previously calculated. These sub- 
routines are called twice by WASP for saturation properties, once with DL and once with 
DV as input for D. 
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SUBROUTINES TO OBTAIN STATE VARIABLES (KS OPTIONS) 

The subroutines used to obtain the state variables are as follows: 

(1) Subroutine DENS (KU, TAU, P, D, DL, DV, KR): This routine solves equa- 
tion (B13) for the density, given TAU and P in units indicated by KU. The region num- 
ber (KR) is returned, and the density is returned in D for KR=2 or KR=3. For KR=1 
the saturation values are returned in DL and DV. If KR=1 for input and either TAU-0 

or P=0 for input, the saturated value is calculated and returned for the variable which 
was input as 0. 

The solution is obtained by ROOT for subcritical pressures and by SOLVE for sat- 
uration or supercritical pressures. Special initial estimates were found necessary for 
convergence near subcritical temperatures with SOLVE and for the interval used by 
ROOT in the region P > P c and 373. 15 K < T < 453. 15 K (100° C < t < 180° C) 

(2) Subroutine PRESS (KU, TAU, D, P, KR): This routine calculates pressure 
(eq. (B13)) as a function of TAU and D in regions KR=2 and KR=3 and as a function of 
TAU only in region KR=1 (using subroutine PSSS). The result, P, is returned in user’s 

units indicated by KU. The correct value of KR is also returned and the calculation is 
direct. 

(3) Subroutine TEMP (KU, P, D, TAU, KR): This routine solves equation (B13) for 
the temperature parameter TAU, given P and D in user’s units specified by KU. In re- 
gions KR=2 and KR=3, SOLVE is used to obtain the solution. In region KR=1, which is 
either input or determined, TAU is a function of P only and is obtained from subroutine 

TSS by solving equation (B31) for TAU. Subroutine TSS also uses SOLVE. The correct 
KR is returned. 

(4) Subroutine TEMPPH (KU, P, H, TAU, D, DL, DV, KR): This routine solves 
equation (B13) by using equation (B20) for the temperature parameter TAU and density D, 
given P and H as input in user’s units indicated by KU. The double iteration is per- 
formed by using ROOT and ROOTX with function TSHF for regions KR=2 and KR=3. In 
region KR=1, the saturation values are determined for DL and DV by DENS and TAU is 
found by function TSS (using SOLVE). KR is also returned. 

(5) Subroutine TEMPPS (KU, P, S, TAU, D, DL, DV, KR): This routine solves 
equations (B13) and vB24) for TAU and D in the same manner as TEMPPH, using P 
and S as input and function TPSF for the double iteration with ROOT and ROOTX. 


SUBROUTINES TO OBTAIN DERIVED THERMODYNAMIC PROPERTIES 

The subroutines used to obtain derived thermodynamic properties assume that the 
variables TAU and D have been input or previously calculated in the user's units. This 
condition is satisfied in subroutine WASP. When KR=1 is input or has been so deter- ' 
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mined, subroutine WASP makes two calls to each routine, once using DL and once DV 
for input D; and the corresponding saturated variable is output [(HL, SL, etc.), (HV, 

SV, etc.)]. 

These subroutines are as follows: 

(1) Subroutine ENTH (KU, TAU, D, H): This routine calculates enthalpy H in 
user's units (KU) by using equation (B20). 

(2) Subroutine ENT (KU, TAU, D, S): This routine calculates entropy S in user's 
units (KU) by using equation (B24). 

(3) Subroutine CPPRL (KU, TAU, D, CP, CV, GAMMA, C): This routine calculates 
the following in user's units indicated by KU: 

(a) Specific heat at constant pressure, CP, eq. (B27) 

(b) Specific heat at constant volume, CV, eq. (B26) 

(c) Specific -heat ratio, GAMMA, eq. (B28) 

(d) Sonic velocity, C, eq. (B30) 

In addition, the first partial derivatives of P are calculated and returned in COMMON/ 
PARTLS/PTV, PDT in the units of KU=1 only. PTV is equation (B16) and PDT is equa- 
tion (B14). 


SUBROUTINES TO OBTAIN TRANSPORT PROPERTIES 

The three routines used to obtain the transport properties assume that the input 
variables for pressure and density and the temperature parameter t are all available 
in user’s units. They are called twice by WASP for saturation conditions, once with DL 
and once with DV as input for density DIN. 

(1) Subroutine VISC (KU, KR, TIN, PIN, DIN, SVISC): This routine uses TIN, PIN, 
and DIN as input in user's units KU. Dynamic viscosity, SVISC, is calculated by using 
one or more of equations (B32) to (B39), depending on the region of the input variables 
as shown in figure 1 and explained in appendix B. All calculations of dynamic viscosity 
are direct evaluations of curve fits. 

(2) Subroutine THERM (KU, KR, TIN, PIN, DIN, EXCESK, TCOND): This routine 
uses TIN, PIN, and DIN in user’s units KU to calculate the thermal conductivity TCOND 
in user's units KU. An optional coding section calculates the critical excess thermal 
conductivity associated with the critical anomaly in the PVT region, 0. 6 < p/p c < 1. 4 
and 0. 9 < T/T c < 1. 1. See also references 12, 13, and 15 and the subroutine listing 
in appendix E. 

The equations used for thermal conductivity are (B40) to (B59) for the different re- 
gions as shown in figure 2. The equation for region IV (et;. (B52)) is iterative. The 
thermal conductivity for the other regions is calculated by direct evaluation of curve fits. 
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(3) Subroutine SURF (KU, KR, UN, SURFT): This routine uses TIN, the input 
emperature parameter, in user's units, to calculate both the surface tension of liquid 
water and the Laplace constant. The calculated surface tension is returned in SURFT 

COMMON/LAPI A°r/r, n r <ALC) ‘ S re ‘ Ur ” ed COMMON /LAPLAC/ALC. The statement 
MMON/LAPLAC/ALC must appear in the user's calling routine if the Laplace con- 
stant is desired. 
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APPENDIX D 


MODULAR DESIGN OF WASP 

A user with limited core storage or with specific property needs may wish to use 
only parts of WASP. The subroutines have been coded so that most of the subroutines 
corresponding to the "KP option” requests may be removed without causing errors in 
logic or calculations. Table V indicates which routines are absolutely necessary and 
which are optional. The conditions for removal must be strictly followed. For sim- 
plicity, the KP options are discussed as though only one option was being requested. In 
reality, the input variable KP is always the summation of the KP option variables. To 
modify a statement number in subroutine WASP, simply replace it with a continue 
statement of the same number. For example, to remove the viscosity option, remove 
subroutine VISC. In subroutine WASP, alteration would read as follows: 

160 CONTINUE 
170 D0175 1=1, 32 

If the user wishes to omit many options, he should rewrite subroutine WASP for effi- 
ciency. 
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APPENDIX E 


PROGRAM LISTING AND FLOW CHART FOR SUBROUTINE WASP 

PROGRAM LISTING 


S1BFTC AQUA 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE WASP I XS * KP ,TT • P *0 , H • KR > 

KEYFS KENNAN HILL MOORE EQUATION OE STATE FOR WATER 

VERSION MARCH 1,1972 

COMPUTE THE STATE RELATIONS AND THERMOD YNAH IC ANO TRANSPORT 
PROPERTIES OF WATER GIVEN TEMPERATURE TT , PRESSURE P, 

DENSITY 0* OR ENTHALPY H * OR ENTROPY S. STATE RELATIONS ARE 
SPECIFIED BY K-S* THERMODYNAMIC ANO TRANSPORT PROPERTIES 
ARE SPECIFIED BY KP . IF KR IS RETURNED OR SPECIFIED AS l* 
PROPERTIES ARE- COMPUTEO AT SATURATION. 


COMMON/ PROPTY/ KU.OL • OV ,HL , HV , S* SL,SV,CV,CVL*CVV,CP ,CPL ,CPV, GAMMA, 
1GAMMAL ,GAMMAV«C,CL ,CVP,MU * MUL , MUV ,K,KL,KV« SI GMA , EXCL , E XC V , E XC E SK 
REAL M U . MOL « MU V • K ,KL , K V 

COMMON /CHECKS /0CH1, DC H2 . PCH 1 * PCH2, PCH3 , TCH1 • TCH2 * TC H3 ,DST , T ST , H 

ISCHlt HSCH2 


C lALi IS TEE TEMPERATURE PARAMETER USED IN THE EQUATION OF STATE 
C TAU IS EQUIVALENT TO T IN THIS SUBROUTINE 

C 
C 

OiFFKSICN K PCI I 32 I, KPC2I32), KPC3 1 32) ,KPC4 1 32 ) 

DATA K PC 1 /2, 3, 6, 7, 10, 11 » 14,15 ,18, 19 ,22 ,23 ,26, 27 ,30 » 3 1 , 34, 35, 38, 
139,42.43*46,47,50,51,54,55 ,58,59,62*63/ 

OATA KPC2 /4, 5, 6, 7, 12, 13, 14, 1 5 , 20, 21 , 22 , 23 ,2 8, 29 *30 , 31 , 36, 37,38, 
139,44.45, 46,47, 52, 5 3,54,55, 60, 61, 62, 63/ 

CAT A K PC 3 /8,9, 10,11,12,13,14, 15 , 2 4, 25 , 26 , 2 7 ,2 8* 29 ,30, 31,40,41 ,42, 
143.44, 45, 46, 55, 56 « 57, 58* 59 ,60,61 ,62, 63/ 

DATA KPC4 /16,17, 18, 19,20,21,22,23,24,25,26,27,28,29,30.31, 
148,49.50, 51 ,52, 53,54,55, 56, 57, 58.59,60,61 ,62 ,63/ 

T*TT 

IF (TT.GT.O.) T * 1000 « /T T 
GC TO I 10,20,30,40,45 1 ,KS 
r 

C COMPUTE DENSITY 

f 

10 CALL DENSIKU.T ,P,0,0L,0V,KRJ 
IF ( TT .EQ. 0.0) TT* 1000 • /T 
GO TO 50 
f 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 
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C COMPUTE PRESSURE 

C 

20 CALL PRES SI KU.T .D.P.KR ) 

GO TO 50 
C 

C COMPUTE TEMPERATURE 

C 

30 CALL TEMPIKU.P.OtT.KR) 

TT-IOOO./T 
GO TO 50 
C 

C COMPUTE TEMPERATURE AND DENSITY GIVEN PRESSURE AND ENTHALPY 

40 CALL TEMPPHlKU.PtHtT.DtDL tDV.KR) 

TT-1000./T 
GO TO 50 

45 CALL TEMPPS I KUt Pt S .T *D* DL « OV • KR > 

TT-1000./T 

50 IE IKR*NE.1.0R. IKS.EO. 1«GR.KS. GT.3 > > GO TO 55 
C 

C OBTAIN SATURATION DENSITIES DL AND DV FOR KS*2 AND KS*3 CALLS WHEN 
C KR* 1 
C 

CALL OENSIKU.T.P.O.OLtOVtll 

55 IF (MO0IKP.2D60, 70.60 
C 

C COMPUTE ENTHALPY 

C 

60 IF IKR.EO.ll GO TO 65 
CALL ENTHIKUtT.O.Hl 
GO TO 70 

65 CALL ENTHI KU.T . OL .HL > 

CALL ENTH IKU.TtDV.HV) 

70 00 80 I* It 32 

IF UP-KPCtf I) U10.100.80 

80 CONTINUE 
GO TO 110 

c 

c COMPUTE ENTROPY 

C 

100 IF I KR «EQ . 1 1 GO TO 105 
CALL ENTIKU.T.D.S1 
GO TO 110 

105 CALL ENT I KUt T t DL .SL ) 

CALL FNTIKU.T.DVtSV* 

110 00 120 1*1.32 

IF! KP— KPC2I ID 14 Ot 130.120 
120 CONTINUE 
GO TO 140 
C 

c COMPUTE SPECIFIC HEATS AND GAMMA AND SONIC VELOCITY 

C 

130 IF I KR .NE • D GO TO 135 

CALL CPPRL I KU.T .DL.CPL .C VL .GAMHAL . CL 1 
CALL CPPRL I KU.T. DV.CP V.C V V .GAMMAV . CVP 1 
GO TO 140 

135 CALL CPPRL I KU« T .D.CP .CV. GAMMA . Cl 
140 on 150 1*1.32 

IF I KP-KPC3I ID 170. 160. 150 
150 CONTINUE 
GO TO 170 
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43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 
70 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

in 

101 

102 

103 
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u u u 


C 

C 


160 


COMPUTE VISCOSITY 


C 

C 

C 


160 IF IKR.NE.n GO TO 165 

CALL VI$C(KU*KR*T +P«DL*MUL > 
CALL V I$C(KU*KR*T #P * D V *MU V ) 
GO TO 170 

165 CAIL VISCIKU*KR*T * P * 0 * MU ) 
170 00 175 1*1.32 

I F ( KP-KPC4I I * > 190* 180*175 
175 CONTINUE 
GO TO 190 


COMPUTE THERMAL CONDUCTIVITY 


180 IF IKR.NE.l) GO TO 220 

CALL THERM (KU *KR .P.T *DL * EXCL* KL) 
CALL THERM ( KU « KR *P 9 T ♦ DV . EXCV ♦ KV ) 
GC TO 190 

220 CALL THERM C KU* KR *P*T* 0. EXCESK ,K) 
190 I F I KP— 32 1 230*240*260 


COMPUTE SURFACE TENSION 

240 CALL SURF I KU. KR * T* S IGMA 1 
230 RETURN 
END 


104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 
129 
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IIBFTC BLGC 

BLOCK CAT A 

C — VERSION MARCH 1*1972 

DIMENSION A It 10).A2< 10),A3<10) *A4< 10) *A5<10J *A6( 10 I »A7 (10 > 
EQUIVALENCE < A 1 <1> * A 1 1 * 1 >)<( A2< 1 ) t A < l *2D • (A3 ( i> *A < 1 *3 ))*< A4( 1 >* A< 
11<4ID<A5< D«A< 1*5) )*(A6<1) • A < 1 * 6) ) * < A7 ( 1 > *A<1*7) ) 

COMMON /CODE / MESSAGI 16 ) 

COMMON /COF / A< 10*71 

COMMON /CRIT/ RHOCRT *PCRT • TCRT 

COMMON /CCNSTS/ T AUC * RHOA * RHCB «TAUA*E *R 

CGMM0N/CH£CKS/DCH1*DCH2* PCH 1 *PCH2* PCH3 *TCHi*TCH2*T CH3* DST • T ST * 
1HSCH1 * FSCH2 
C 

COMMON /CO SAT /CP SI * CP$2.CPS3.CPS4*CPS5.CPS6. CPS7 
COMMON /StCOF/ S IC1* SIC2*SIC3* SIC4 *SIC5 
COMMON /CON VI /OCONV <51 
COMMON /C0NV2/T CONV < 5 ) 

COMMON /C0NV3/PC0NV < 5) 

COMMON /CON V4/ SCON V ( 5 1 
COMMCN/CON V5/CC0NV< 5 ) 

COMMON/CONV6/HCON V< 5 ) 

COMMON /CGNV 7/ MCONV < 5 > 

COMMON/CON V8/KC0NV< 5) 

COMMCN/CON V9/STC0NVI 51 
REAL MCONV • KCGNV 

OAT A A 1/29 *492937 *-132.139 17*274.64632 *-360. 93828* 342. 1843 1* 

I -244.50042* 155.18535* 5. <1728487* •410.30848* -416.05860 / 

OATA A2/ -5.1985860* 7.7779182* -33.301902* -16.254622 *-17 7.31074* 
1127.48742* 137.46153* 155.97836* 337.31180* -209.88866 / 

OATA A 3/ 6.8335354* -26.149751* 65.326396* -26.181978* 4*0.* 
1-137.46618* -733.96648 / 

OATA A4/ -.1564104* -*72544il08* -9.2734269* 4.312584* 4*0. • 

1 6.7874583* 10.401717 / 

OATA A 5/ -6.3972405* 26.409282* -47.740374* 56.32313, 4*0.* 

1 136.87317* 645.8188 / 

OATA A 6/ -3.9661401, 15.453061* -29.14247* 25.568796* 4*0.* 

1 79.84797* 399.1757 / 

OATA A 7/ -.69048554* 2.7401416* -5*102807* 3*9636085* 4*0.* 

1 13.041253* 71.531353 / 

OATA CPS1.CPS2*CPS3*CPS4*CPS5*CPS6*CPS7/ 

1 0*2930437 OF >01 * - 0.23095789E+04* 0. 34522497E-01 * 

2 —.1362 1289E-03* 0, 25878044E-06 • -.24709162E-09 * 

3 0. 9593764 6E- 13 / 

OATA 0CH1 * DCH2 • PC Hi * PC H2 • PCH3 * TCH1 * TCH2 * TCH3 *0 ST *T ST* 

1HSCH 1* HSCH2/0. * 1 .04539 « .00 06 1 1 *22.089*100. •• 55555556* 1.544912* 

1 3.660900 *. 8*400**0. *4700* / 

DATA MESSAG / 96H THER MOOYNAM 1 C AND TRANSPORT PROPERTIES FOR MATER 
1PC*218.07ATM*TC* 647.29 K*R0C=.3L7 G/CC / 

DATA RHCCRT*PCRT*TCRT/.317. 22*089* 647.286 / 

OATA SIC1*SIC2*SIC3.S1C4*S IC 5/ 1855.3865* 3.2 78642* 

1-0.00037903* 46*174* -1.02117 / 

OATA T AUC* RHOA* RHQ8* TAUA* E »R/ 1.544912* .634* 1«* 2.5* 4.8*.4615l/ 
OATA ICON V /1..1. ..55555555 .2*1.0 / 

OATA PCCNV /l. *9.8 69 2 32 7*145. 03824 3*2*1*/ 

OATA 0CCNV/2*1 * . 62*4283* 2* l * ✓ 

OATA SCCNV/2*1.*0. 238849*2*1./ 

OATA CC0NV/2*l«»0.0328l*2*l./ 

OATA HCCNV/2*1 • *0.429929*2*1. / 

OATA MCONV/ 2*1 .0. .671 ‘669 9E-1 * 2*1 . 0/ 

OATA KCCNV / 2*1 *E— 2* 1.6C604 4E— <* *2*l.fc-2/ 

OATA STCON V / 2*1 . * 6.852 1 766E— 5 • 2*1 ./ 

COM MO N/XMINDS/XM1 (7) 

OATA XM1/1.*2.« B.*4..5.*6.*7./ 

END 
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: ROOT 1 

771177 — VERSION 

SUBROUTINE ROOT* C XO* X2 .FOFX * FUNC* XI) 


SOLVE FOR XI SUCH THAT FUNC I XI) 
BETWEEN XO AND X2 


FOFX, WHERE XI LIES 


COMMON /CHECK 2/KOUNT 

TOL-I.E-5 

XXO * XO 

XX2 ■ X2 

FO » FUNC I XXO) 

F2 * FUNCIXX2) 

A* I EOF X-FO ) / ( F 2— FO) 

IF I A ) 1007*120* 120 
120 IF IA-1.) 130*130*1008 
130 IF 1 FOFX— 0 • ) 80*70*80 
70 ASSI6N 100 TO JUMP 
GO TO 90 

80 ASSIGN 110 TO JUMP 
90 X « I X X0«XX2 ) /2 • 

MOUNT * 0 
150 XI - X 

MOUNT « MOUNT ♦ 1 
A * FOFX - F2 
FX - FUNC I X ) 

FXL»F0*1X— XXO)* IF2— FO) /( XX2— XXO) 

B«ABSI I FX— FXL 1 / IF 2— FO) ) 

IF I A* ( FX— FOFX ) *LT • 0.) GO TO 1001 

XXO » X 

FO-FX 

IF IB-.3I 10*20,20 
2C X « (X+XX2 )/2* 

GO TO AO 
1001 XX2 - X 
F2 » FX 

IF IB-.3) 10.30,30 
30 X a lXX0*X)/2. 

GO TO AO 

10 X~XX0+(F0FX-F0)*( XX2-XX0) /(F2-F0) 

AO IF I ABSI I X— X 1 ) / X 1 -TOL ) 50,1000,1000 
50 GO TO JUMP. 1 100,1 10) 

100 IF IABSIFUNCIX) 1-TOLAlO. >60, 1000, 1000 
110 IF IABS 1 1 FOFX— FUNCI X ) ) /FOFX )— TOL ) 60,1000.1000 
1000 IF I MOUNT «GT *A0 ) T0L«TCL*1C. «>u, lUUO. 1000 

IF I MOUNT *GT *60 ) TOL«TOL*10. 

IF (MOUNT » GT *80 ) T0L*T0L*10« 

IF (MOUNT*! T« 100) GO TO 150 
160 WRITE 16.170) XI, X 

6C X1*X 
RETURN 

1007 XI * XO ! 

GC TO 1 AO 
10C8 XI = X2 
l AO WR ITEI6.1A1 ) 

1^1 FORMAT ( 1H0 . 2AH SOLUTION OUT OF RANGE 
RETURN 
ENC 
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IIRF1C 8G0T2 

SUBROUTINE ROOT 


4 XO.X2 .FOFX .FUNC. XI) 


f 

C 

c 

c 

c 

f 


VERSION 2/1/72 

SAME AS ROOTX - NEEDED TO PREVENT RECURSION 
SOLVE FOR XI SUCH THAT iFUNCCXil - FOFX. WHERE XI LIES 
BETWEEN XO AND X2 


COMMON /CHECK I /MOUNT 

TOL-l.E-5 

XXO * XO 

XX2 * X2 

FO * FUKC ( XXO ) 

F2 * FUNC ( XX2 > 

A«IFCFX-F01/tF2-F01 
IF ( A I 1007 • 120 • 120 
120 IF (A-l.l 130.130*1008 
130 IF 4F0FX-0.I 80.70.80 
70 ASS IGN 100 TO JUMP 
GO TO 90 

80 ASSIGN 110 TO JUMP 
90 X * (XX0 + XX2 1/2. 

KGUNT * 0 
150 XI * X 

MOUNT * MOUNT ♦ 1 
A * FOFX - F2 
FX - FUNC I X I 

FXL*F0*4X— XXO I* 4F2-F01 /I XX2-XX01 

B*ABS4 4FX-FXL1/4F2-F011 

IF I A* I FX-FOFX I .LT. 0.1 00 TO 1001 


20 

1001 

30 

1C 

40 

50 

1O0 

110 

1000 


16C 

170 


XXO - X 
FO-FX 

IF ( 0— .3 1 10.20.20 
X * I X+XX2 1/2. 

60 TO 40 
XX7 * X 
F2 * FX 

IF 1 8- .3 1 10.30.30 
X » IXX0+X1/2. 

GO TO 40 

X«XX0*(FCFX-F0I*4XX2-XX01/4F2-F01 
IF 4A8SI4 X-Xl 1/XI-TOL 1 50.1000.1000 

TO JUMP. 4 100.1101 

(ABS I FUNC ( X 1 1-TCL410. 160*1000. 1000 
( ABS ( 4 FOFX— FUNC ( X 1 1 /F0FX1-T0L 1 60.1000. 1000 
(MOUNT .GT.401 T0L»TCL*1G. 

( MOUNT. 6T .601 IGL*TOL*10. 

4M0UM.GT.801 T0L«T0L*10. 

(MOUKT.LT. 1001 GO TO 150 

FORMAT ^ t lHi?79HAN*iTERAT ION HAS BEEN TERMINATED AT 100 
1 TFE LAST TWO VALUES WERE .3G15.51 


GO 

IF 

IF 

IF 

IF 

IF 

IF 


6C X 1 =X 
RETURN 

1003 XI = XO 

GO TO 140 
1008 XI * X2 

140 WR ITE46.1411 

141 FORMAT ( 1H0.24H SOLUTION OUT OF RANGE 
RETURN 

ENC 


ITERATIONS. 


1 

2 

3 
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17 
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21 
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27 

28 
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36 
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30 

39 
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48 
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SIBPIC SC L V 

c VERSION 2/1/72 — 

FUNCTION SOLVE I X I #F# OF I 

C 

C NEwTCN-RAPHSON I TER AT (OK GIVEN AN INITIAL ESTIMATE XI 

C AND THE FUNCTIONS F ANO OF 

CGPMGN /CHECK 1 /NI 
TOl-l.F-5 
N 1*0 
XO*XI 
XN-XI 
10 XOO*XO 
XG*XN 

XN*XO— F I XO 1 /OF ( XO I 
NI«NI+1 

IF IABSI IXN-XOI/XNI-TOL 1 70#20*20 
20 IF IN I #GT *40 ) TOL*TOLMO. 

IF INI *GT*60) TCL~TOL*iO. 

IF INI *GT *601 TCL*TOL*lO* 

IF INI— 100 ) 30* 50*50 
JO IF I ABSI 1 XN-XOG l/XNl— TOL 1 40*10*10 
4C XN* I XO+XN 1/2* 

GC TO 10 

50 WRITE I 6*60 > XOO*XO*XN 

60 FORMAT <1HL*B1HAN ITERATION HAS BEEN TERMINATED AT 100 ITERATIONS* 
1 THE LAST THREE VALUES WERE *3G15*5> 

70 SGLVE-XN 
RETURN 
ENO 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
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22 

23 

24 
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SlBFTC scfec 

FUNCT ION CHECK (KluKK.T) 


( 

C 

c 


VERSION MARCH 1,1972- 


CCMMCN/CCNV 1/0C0NVI 5) 

CfiMMGN/C0NV2/TC0N VI 51 
COMMON /CON V3/PC0NVI 51 
COMMON/ IERROR/ IROUT 

COMMON /CHE CKS/DCH l «DCH2« PCHi « PCH2, PCH3*TCH1 » TC H2 , TCH3* OST »TST » 


1HSCH 1 * FSCH2 

0 IMEK S ION FM1<9>* FM2I91* FM319), FMTC 9) • RQUTClil 

OA1A FM1 /51HI1H ,61 2*4 »31HIS OUT OF RANGE FOR T IN SUB*-,A6 > 

^OATA FM2 /51HI1H »G12*4 *3iHIS OUT OF RANGE FOR P IN SUB.-*A6 I 

*DATA FM3 /51HUH .G12.4 •31HIS OUT OF RANGE FOR 0 IN SUB.-.A6 J 

DATA ROUT /4HDENS, 5HPRESS ,4HTEMP ,4HENTH , 3HENT , 6HTEMPPH, 6HTEMPPS 
1 ,5HCPPRL*4HV I SC ,5HTHERM* 4HSURF / 


C CONVERT TEMPERATURE T TO DEGREE $ KELVIN AND CHECK 

C FOR OUT OF RANGE# UNITS ARE SPECIFIED BY KU# IF KR 

C IS SPECIFIED AS l. T IS CHECKED FOR OUT OF SATURATION 

C RANGE# 

C 

ENTRY TCHECK ( KU* KR , T I 
CHECK* 1000 #*TCON V ( KU > / T 
CHI* 1000 • / TCH3 
CH2M000./TCH2 
CH3* 1000* / TCH1 
KCCE-l 
OC 1 J* 1 1 9 
1 FMT I J l*f M 1 1 J > 

GC TO 10 


CONVERT PRESSURE TO MN/M**2 ANC CHECK 
FOR OUT OF RANGE# UNITS ARE SPECIFIED BY KU# IF KR IS 
SPECIFIED AS 1, P IS CHECKED FOR OUT OF SATURATION 
ENTRY PCHECK f KU* KR * P 1 
CHECK*P/PCONV(KU> 

CHI* PCHI 
CH2* PCH2 
CH3* PCH3 
K0DE*0 
00 2 J* 1*9 
2 FMT 4 J ) *FM2 C J I 
GO TO 10 


CONVERT DENSITY TO G/OC ANC CHECK 
FOR CUT OF RANGE# UMTS ARE SPECIFIED BY KU. 

ENTRY CCHECK ( KU *D ) 

CHECK *0/DCCNV I KU > 

CHl*DCHi 
CH3*CCH2 
KCCE*0 
DC 3 J* 1 • 9 
3 FMTI J1*FM3I J> 

GO TO 20 

10 IFIKR.FG.n GO TO 30 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
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23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
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20 I E (CHECK *L T *CH1 1 GO TO 40 
I f (Cf ECK.GT «CH3 1 GO TO 40 
25 If (KOCE • EQ . 1 1 CHECK*T/TC0NV(KU1 
RETURN 

BO If ( CHECK. L T .CHI > GO TO 40 
IF(CHECK.LE.CH21 GO TO 25 
40 WRITE! 6*FMT I CHECK* ROOT II RCtT > 

GO TO 25 
ENC 


61 

62 

63 

64 

65 

66 

67 

68 
69 


AIBFIC Sit B 1 

SUBROUTINE CRUST (01 

C VERSION MARCH 1*1972* 

COMMON /CCNS TS/ T AUC *RHOA *RHOB « TAUA • E *R 

COMMON /QAUX / RBOIf 4 6) *RAOIf I 8) *ER*EO» TAOlf I 7) 
COMMON /COf/ A I 10*71 
COMMON /OS 1 /SUM 1(7) 

R AC If (II* | # o 
RACIf (21 * 0- RHOA 
RBOIf (11 * 1.0 
RBCIf (21 * 0- RHOB 
00 1 I* 3.8 

RBCIf (11 * RBOIf 11-11# RBO I f ( 2 1 
1 RAO If (II* RADIf (1-11 ARAD I f ( 2 1 
EO * FAD 
f R* 1.0/ EXf(ED) 

SUM I ( 1 1*0.0 
00 4 1*1*8 

4 SUM 111 I* SUM 1(11 ♦A (I* 11 ARAD If ( II 

SUMI ( 1 1*SUMI( 1 1 *ERA( A< 9* 1 1 *A ( 10*11A 0 1 
00 6 J*2* 7 
SUM 1 1J 1*0.0 
00 5 1*1*6 

5 SUM 1 1 J 1*SUM I( J1 *A ( I • J1 ARBOIf ( I 1 

SU M I ( J 1* SUM ((Jl+ERA(A(9tJl+AU0*J)ADl 

6 CONTINUE 
RETURN 

ENTRY 0MUST2I T AU1 
TAC If (I 1 * 0.0 
TACIf (2) * 1.0 
TACIf (31 * TA l- TAUA 
00 2 I* 4*7 

2 TACIf (II* TAO I f ( 1-1 |A TAD I f ( 31 
RETURN 
ENC 


1 

2 

3 

4 

5 

6 
7 
6 
9 

10 

11 

12 

13 

14 

15 

16 
17 
16 

19 
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22 
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25 
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4IBFTC 5082 

FUNCTION QCALC1TAU1 

C 

C VERSION MARCH 1*1972 

C THE FUNCTION QIRH0.TAU1 

C 

COMMON/ CHE CKS/OCH 1.DCH2* PCH l* PCH2* PCH3 • TCH1 « TCH2 » T CH3* OST *TST, 
1HSCH1 • HSCH2 

CCMMCN /OAUX/ RBDIFI81* RAO IF! 8> .ER.EO, TA0IFI7) 

COMMON /COF / A( 10* 1 1 
CCMMCN/CS1/SUMI 171 
TSUM * 0*0 
00 4 J*2» 7 

4 T SUM-T SUM*T AO I F ( J »*SUM I ( J > 

OCALC-SUMIt !>♦( TAU-TCH21*TSUM 

RETURN 

ENC 


1 

2 

3 

4 

5 

6 
7 

a 

9 

10 

u 

12 

13 

14 

15 

16 


IIBFTC S083 

FUNCTION OOTAITAUI 1 

C VERSION MARCH 1*1972 2 

C 3 

C PARTIAL OER OF 0 — - PQ/PRHC 4 

COMMON /O AUX/ R80 IF! 8> *RAO IFt 8 » , ER .60. TAOI F I 7» 5 

COMMON /CO F/ AI10.7) 4 

COMMON /CONSTS/ T AUC .RHOA.RHOB , TAUA.E ,R 7 

COMMON / XM INUS/XM1 17) 8 

COMMON /OS 2 /SUM 1(7) 9 

EQUIVALENCE I SUM 1 1 1 1 • SUM) 10 

1 TSUM-0.0 11 

00 2 J*2» 7 12 

2 T SUM-T SUM 4 T AD I F I J >*SUM II J 1 13 

0DTA-SUM4I TAU-TAUC1*TSUM 14 

RETURN 15 

ENTRY OCTIO.TAUI 14 

SUM-0. 0 IT 

OC 10 1-2.8 18 

10 SUM-SUM-tXM l(I-ll*A(I.l 1*RAC IF I I— 1 A 19 

SUM-SUM*ER«IA( 10* il-E* t A ( 9, 1 1+AI 10,11*01 1 20 

OC 15 J-2.7 21 

SUM 1 1 J 1-0.0 22 

OC 12 1-2*8 23 

12 SUM I i J 1-SUM IIJl + XMH l—ll*All.J )*RBDI FI I— l) 24 

SUNK 01— SUM 1 1 J1*ER*I A( 10* Jl-6* I A 19 * J1*AI 10. J >*D) 1 25 

15 CCNT INUE 26 

GO TO I 27 

FNC 28 
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M8FTC St 84 
^ FUNCTION 

C 

C 


OTCI TAU1 


•VERSION MARCH I • 1972- 


18 


PARTIAL DER OF 0 — — PQ/PTAU 

£2S2!/£“!Em , " H0 *’ 1 RHCa - •' E -« 

COMMON /XM INLS/XM147* 

T SUM1 * 0*0 
TSUM2 * 0*0 
OC 18 J* 3 f 7 

TSLM1«TSUM1+XM1<J-2)*TA0IF< J-l I*SUMI(J1 
TSUM2~TSUM2*TA0JF( J I ♦SUM I ( J1 
T SUM2*T SUM2+SUM I ( 2 ) 

OTO*T SUM2 + f TAU-TAUC J*T SUM1 

RETURN 

ENC 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 


S 16 FTC SUBS 

FUNCTION C212DITAUJ 

~ • — VERSION MARCH 1 * 1972 

-PAR T I AL CER OF 0 P20/PTAU2 

cK; ; c “oT; :r;i:i; , -** o,f,,, - E "- 6D - i “ iF, » 

CCMKON /COASTS/ T AUC ♦ ft hO A , ft HC£» , TAUA, E ,R 
COMMON /OS 1 /SUM 1(7) 

CCftPON/XH INUS/XMl 4 7> 

TS0M1 * 0.0 
IS0M2 « 0.0 
or 2 J*3 .1 

TStMlaTSUHI*XMl( J-2)*TAOIF4 J-i )*SUM( ( J) 

IF IJ.EC.3) GO TO 2 

2 CChllNUE liHi * XM1U “ 2> * XW1< J - 3, * TA0 « F < J-2J*SUHIIJ) 

Q21 2C*2 . 0*T SUM 1 +( TAU-IAUC )*TSUM2 

RETURN 

ENC 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
19 
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IIBFTC SUB6 

FUNCTION 02DTA4TAU1 


C 

C 

c 


VERSION MARCH 1*1972 

• PARTIAL OER OF 0 P20/PRH0— PTAU 

COR MON /QAUX/ RBO I F 4 8 1 *R AO IF 4 8) • ER.ED tTAOlF (7 1 
COMMON /COF/ A I 10*7) 

COMMON /CONSTS/ T AUC •RK3A .RHCB * TAUA . E *R 
COMMON /OS 3 /SUM 4 6 1 
COM MON /XM INUS/XMI 47) 

1 T SUM 1*0*0 
TSUM2*0 *0 
00 10 J*3. 7 

TSUM1*TSUM1+XM1 4 J-21*TADIF 4 J- i l^SUM 4 J-l 1 
10 TSUM2*TSUM2*TADIF< J1*SUMI J- 11 
T SUM2*T SUM 2 ♦SUM 4 1 1 
TSLM1*TSUM1*(TA4J-TAUC1 
02CTA*TSUM1*TSUM2 
RETURN 

ENTRY 020 T I 0* TAU) 

00 20 J*2* 7 
SUM 4 J- 1 1*0 • 

00 15 1*2*6 

15 SUM! J-ll*SUM4 J-11>XM1( I- 1 ) *>A 4 1 • J ) *RBD I F 4 I-l> 

SUM4 J-l 1*SUMI J-1 1*ER*< A C10»J)-E*4A49*JM-A410*J)*D)1 
20 CONTINUE 
GC TO 1 
ENO 


1 

2 

3 

4 

5 

6 
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a 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
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20 
21 
22 

23 

24 

25 
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6I8F1C SUB7 

FUNCTION 02D2TA4TAU) 


C 

C 

c 


VERSION MARCH 1*1972- 

PARTIAL OER OF 0 P20/PRH02 

COMMON /OAUX / R8D IF 4 8 ) *RAO IF 4 6 1 *ER «ED* TAOIF ( 7> 

COMMON /COF/ A 4 10*71 

COMMON /CONSTS/ TALC* RHOA* RHG8# TAUA# E* R 
COMMON /0S4 /SUM 1471 
CO MMON/XM INUS/XMI 4 7 ) 

1 TSUM*0.0 
00 5 J*2.7 

5 T SUM*T SUM* TAOIF 4 J 1*SUM 1 4 J 1 

02C2TA*SUM 1 4 1 1 ♦ 4 T AU-TAUC1 *T SUM 
RETURN 

ENTRY 0202 T 4 0* TAU 1 
SUM III 1*0.0 
00 3 1*3.8 

3 SUM II 1 1*SUM 1 1 1 1+XM1I 1-1 1*XM14 1-21 ♦A 4 I * 1 )*RAOIF 1 1-21 
SUM III 1*SUMII 1 1 ♦ER*4-E*A4 10 . 1 ) * 4 2. O-ED 1 *E*E* A4 9* 1 1 1 
OC 10 J*2* 7 
SUM 1 1 J 1*0*0 
00 8 1*3.8 

8 SUM 1 4 J I* SUM IIJI+XM1! 1- 1 1*XM14 1-2 1 * A 4 I * J 1*RB0 IF 4 1-21 
SUM 1 1 J 1* SUM 1 1 J I +ER*I -E* A 4 10 • J 1* ( 2 •O-EDl *A (9 * J > 1 
10 CONTINUE 
GC TO 1 
ENO 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 
27 


50 



9° o nnnno 


IIBFTC SUBS 


SUBROUTINE PSSSIPSSI 

VERSION MARCH 1,1972 

COMPUTE SATURATION PRESSURE PSSS IN BARS AS A FUNCTION OF T IN DEGREES 
C AND RETURN ANSWER IN PSS IN MN/M**2 

COMMON/COSAT/ CPS1 ,CPS2,CPS3,CPS4,CPS5, CPS6,CPS7 

THE T IN THE COMMON BEN017 IS REALLY TAU 

COMMON/TPARAM/T 
DIMENSION CTIPSI6) 

DATA CTIPS / .31602383E-03 , 1.00044775, -0.46487771E-05, 

1 0 .69431 852E-08, 0. 15621 197E-12 , 1.00043357 / 

TSC » 1000. /T -273.15 

— CONVERT TSC (THERMODYNAMIC CELSIUS TO INT. PRACTICAL SCALE (Cl WHICH 

IS USEO IN SATURATION EQUATION 

IF (TSC «GE. 9.996) €0 TO 9 
TS ■ CTIPS(6 1 * TSC 
GO TO 10 

9 TS « ( ( (CTIPS(5)*TSC«CTIPS(4) )*TSC ♦ CTIPS( 3) )*TSC ♦CTIPS(21) 

1 *TSC ♦ CTIPS(l) 

10 TS-TS+273.15 

PSS»10.**( ( ( ( (CPS7*TS4CPS6)*TS*CPS5)*TS+CPS4)*TS+CPS3)*TS*CPS2/TS+ 
1CPS1 1 

PSS-PSS/10.0 

RETURN 

ENO 


1 

2 

3 

4 

5 

6 

7 
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10 
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$ 1 6FTC STJE9 

FUNCTION TSS I PS I I 

C - VERSION MARCH 1,1972 2 

C 3 

C COMPUTE SATURATION TEMPERATURE IN DEG C AS A FUNCTION OF PRESSURE A 

C IN BARS ANO RETURN ANSNER TSS AS TAU IN KELVIN**- 1 5 

C 6 

C0NN0N/CHECKS/0CHim,0CH2,l»CHl,PCH2,PCH3,TCMl,TCH2,TCH3,0ST,TST,H 7 

1SCH1.HSCH2 8 

CGMNON/COSAT/ CPS1 ,CPS2,CPS3 r CPS4,CPS5,CPS6,CPS7 9 

C0NM0N/BEN09/A1 , A2, A3, A4,A5 10 

OIMENSICN CTTI6) 11 

OATA CTT /-.30733645E-03, 0 .99955209, 0.46490458E-05,-.69336443E-08 12 

1,-0. 18086305E- 12 , .99956701 / 13 

EXTERNAL TSSF.OTSSF 14 

PSI*PS*10.0 15 

Al«CPSl-AL0G10( PS11 16 

A2»5.*CPS7 17 

A3«4.*CPS6 18 

A4»3.*CPS5 19 

A5»2.*CPS4 20 

TESTM *( 1000./TCH2 -20.0 ) 21 

T SS» SOLVE ( TESTM, TSSF, CTSSF ) 22 

TSS»TS$-273.15 23 

C CONVERT THE CALCULATED SATURATION TEMP. FROM INT. PRACTCAL SCALE 24 

C (Cl TO THERMOOYNAMIC CELSIUS SCALE 25 

TSIP-TSS 26 

IF ITSS .GT. 10.) GO TO 9 27 

TSS » CTT 161* TSS 28 

GO TO 10 29 

9 TSS * IUCTT<5)*TSS ♦ CTT(4))*TSS ♦ CTT(3))*TSS ♦ CTT(2))*TSS ♦ 30 

1 CTT ( 1 ) 31 

10 TSS » 1000./ (TSS+273. 15 ) 32 

RETURN 33 

ENO 34 


AIBFTC SUB 10 

FUNCTION TSSF(TSS) 1 

C VERSION MARCH 1.1972 2 

C 3 

C FUNCTION USED TO SOLVE FOR SATURATION TEMPERATURE TSS 4 

0 GIVEN PRESSURE 5 

C 6 

CONMON/COSAT / CPS 1 .CPS2 .CPS3.CP S4.CPS5.CP S6 .CPS 7 7 

CCMM0N/8EN09/A 1.A2.A3. A4* AS 8 

TSSF*< I <ICP$7*ISS*CPS*>*TSS*CPS5)*TSS+CPS4>*1SS*CPSJ)*TSS*CPS2/ 9 

IT SS+A 1 10 

RETURN 11 

ENTRY OTSSFITSSI 12 

C 13 

C DERIVATIVE OF FUNCTION LSED TO SOLVE FOR SATURATION 14 

C TEMPERATURE TSS GIVEN PRESSURE 15 

f 16 

TSSF*I I ( A2*TSS-»A31*TSS+A4 J*TSS+AS >*TSS*CPS3-CPS2/ (TSS + TSS ) 17 

RETURN 18 

ENC 19 
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1 18 FTC PRES 1 

SUBRCUTINE PRESS(KU,T.O,P.KR) 

VERSION MARCH 


C 

c 

c 


COMPUTE PRESSURE P GIVEN TEMPERATURE T ANO DENSITY n 

S'" *" « u - ««■ «« » E fS.«S o" °* 

fuwJmJ OF I CC " PUIEI1 *' S*'U«»TIC, «s » 

COMMON /CONV3/PCONVIS) 

COMMCN/TPARAM/TS 

COMMON /CCNSTS/ TAUC,RHOA,RHOB,TAUA*E ,R 
1HSC»u!‘eSCH2 KS/0CW 1 ,0CH2 ’ PCH 1 ’ PCH2 ’ PCH3 * TCH1 ’ TCH2 * TCM3 » 0ST » T ST , 

CCMRGN/IERRGR/IROOT 

IRGUT«2 

T$*TCHECMKU*KR*T J 
OETERHINE REGION 


IF UR-11 10* 70* 10 
10 OS*OCfiECK ( KU* D I 

IF I KR *GT. 1) GO TO 60 
IF ( T S-TCH2 1 50*50*20 
20 CALL OENS(l*T$*ZE*ZE*OSL f DSV*il 
IF ICS-OSLI 30*60*40 
30 IE IGS-GSVJ 50 * 6C *60 
40 KR«2 

GO TO 60 
50 KR* 3 

GO TO 60 
C 

C REGION 1 


60 KR*I 

70 CALL PSSS(PS) 

GO TO 90 
C 

f REGIONS 2 ANO 3 

C 


8C CALL ONUS T ( C S) 

CALL 0*UST2<TS> 

FS* 1000 **R AOS/ TS4( 1 • ♦OS* lOCALC(TS) 
90 P«PS*PCGNVUU> 

RETURN 

ENC 


+GS*QOT(DS,TS>)> 


10 

U 

12 

13 

14 

15 

16 
17 
16 

19 

20 
21 
22 

23 

24 

25 

26 
27 
26 

29 

30 

31 

32 

33 

34 

35 

36 

37 
36 

39 

40 

41 

42 

43 

44 

45 
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IIBFT C CEN5 1 

SUBROUTINE CENSCKU*T*P*0*OL*DV*KR> 

c - VERSION MARCH 1.1972 

C COMPUTE OEMS IT V D GIVEN TEMPERATURE T ANO PRESSURE P. 

C UMTS ARE SPECIFIED BY KU. IF KR IS RETURNED OR 

0 SPEC I F I EO AS I* THE SATURATED LIQUID ANO VAPOR DENSITIES* 

C DL ANO DV RESPECTIVELY* ARE COMPUTED AS A FUNCTION 

C OF T OR P. THE OTHER VALUE MUST BE INPUT AS 0*0 • 

C 

COMMON /CHECK 1 / N I 
COMMON /CONVI/OCONVi 51 
COMMON /CON V2/TCCN VC 51 
C0MMCN/CCNV3 /PC0NVC5) 

COMMON/ I ERROR/ K ROUT 

COMMON /CRIT/ RHOCRT *PC8T * TCRT 

COMMON/ PS I CON/ S IC1 • S IC2 * SIC3 • S IC4* S ICS 

COMMON /CONSTS/ T AUC *RNOA *RHG8 « TAUA * E *R 

C GHMON / CHECKS / OC H 1*0CH2* PCH1 * PCH2* PCH3 * TCHl * TC H2 * TCH3* DST *T ST * 
IHSCH1 • NSCH2 
COMMON /TP AR AM/ T S 
COMMON /PR HOT/PS* OS* TT 
EXTERNAL CSF* DOSE 
IRGUT* 1 

IF ( KR • EO • I 1 GO TO 70 
TS*TCHECK(KU*KR*TI 
TT * TS 

CALL 0MUST2ITS1 
GO TO 5 

70 IF I T *6T .0.0 1 GO TO 75 
PS* PC HECKCKU*KR*P » 

T S*TSS I PS I 
TT * TS 

IF CT-LE.O.I T® TS *TC0NVCKL1 
CALL GMLST2ITS1 
GO TO 5 

75 TS * TCRECK (KU*KR«T1 
TT * TS 

CALL GMUST 2 ( TS 1 
CALL PSSSCPS) 

IF CP.LE.0.1 P*PS APCONVCKUI 
C 

c DETERMINE REGION 

C 

5 IF ( KR- 1 1 10*80*10 
10 PS* PC HECK < KU* KR * P ) 

IF cps-pch; 11 10*1 10* 100 

100 IF (TS-TCH2U30 * 130* 120 

120 KR*2 

FST1 * 1*0455 

EST2 * RHOCRT 

TEST * 1000* /TS - 273.15 

IF t TFST .GT. 100.1 EST1 * 1.0107054 

IF ( TEST .GT. 180.01 GO TO 121 

IF (TEST #L T .40 • 1 TEST*40. 

ES T2= ( IE STM TEST* I TEST*! TEST* . 1247 671 1 E-09-. 52277795E-07 1 
1 ♦ • 54 79057 1 E-05 1— .696 1 7325E— 03 )♦ 1 • 02202 77 1 

121 CALL RCCT ( EST i * E ST2 • 0 . * C SF * D S 1 
GC TO 150 

130 KR * 3 

FST*RH0CRT*3. 
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CALt ROOT! EST *DCH|*0* *C$F *CS) 

GO TO 150 

110 IF 4TS-TCH2 I 50*50*20 
20 CAIL PSSSIPSS) 

IF I AS SI 4 PSS~PS)/PSS )*1 *E-4) 6C*30,30 
30 IF 4PS-PSS) 50*60*40 

40 KR*2 

41 05**71654566 

IF 4 TS *GT* 1*74474391 OS * *622368 
IFITS *61* 1*94674791 OS * *6474576 
IF 4 TS *GT* 2*02776051 OS * *907441 
IF 4 TS *6T. 2*30666901 0S**961538 
IF 4 TS *GT« 2*6798674) OS » 1*001001 
IF 4 Kfl *£0* 1) GO TO 81 
GO TO 90 
50 KR*3 

OS * PS4TS/I 1000*48) 

GO TO 90 


REG ICN 1 


60 KR* 1 

80 CONTINUE 
GO TO 41 

61 CONTINUE 

OSL*SOL VE4 OSvOSF* OOSF ) 

OS « PS*TS / 4 1000* *R ) 

IF 4 TCH2/TS *GT • *965 IDS** 654RHQCRT 
IF 4 TCH2/TS*GT • *995 ) DS»* 75ARHOCRT 
IF 4 TCH2/T S *GT • *999 ) DS»* 65*RHCCRT 
IF 4 TCH2/TS*GT* *9995 ) OS* *90*RH0CRT 
OSV*SOL VE4 OS*OSF* OOSF ) 
0L*DSl*CC0NV4KU ) 

OV*0$V*0CGNVUU> 

RETURN 

C 

C REGIONS 2 ANO 3 

C 

90 OS*SGLVEIOS*OSF*OOSF) 

150 0*OS«OCGNV4KU) 

RETURN 

ENC 
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< [BMC OSFI 

FUNCTION OSFIO » 

C FUNCTION USEO TO SOLVE FOR OENSITV 0 GIVEN TEMPERATURE 

C ANO PRESSURE 

COMMON /CONSTS/ TAUC»RROAtRHOB«TAUA»E »R 
COMMON /PRMCT / PS»OS«TS 

PST ATE* 1000 *RR*0 /T$*l !•♦!) ♦ COCALC ITS) *0 ♦OOTID »TS>II 
OSF-PSTATE-PS 
RETURN 
C 

ENTRY CCSFIC » 

OOSF»iOOO«iS/T S A( 1^+D* 12«0*0CALC CT SI *4 #0*0*0 OT ID tT SI *0*0*02 D2T COt 

1TS1II 

DSf*OOSf 

RETURN 

ENO 
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AIBFTC TENPi 


SUBROUTINE TENP(KU,P, 0,T,KR ) 

VERSION NARCH 1,1972 


S?TV UT IF T KR“lS 0 fp^Jf2?l5 A I« R f ii USERS UNITS GIVEN ASSURE ANO OEN- 
5ITV. IF KR IS SPECIFIED AS l TAU WILL BE A FUNCTION OF PRESSURE ONLY 


C 

C 

C 


C 

C 

C 


C 

C 

C 


C 

C 

C 


COPPON /C0NV2/TC0NVC 5 ) 
COPPON/CHECKS/OCH1 ,DCH2, 
1HSCH1, HSCH2 
COPPON/IERROR/ IROUT 
COPPON /PRHOTT/ PStDStTS 
EXTERNAL TSF*0TSF 
IROUT-3 


PCH1 * PCH2»PCH3 f TCH1, TCH2,TCH3,DST,TST» 


PS«PCHECK(KU,KR,P) 


DETERPINE REGION 


IF IKR-l) 10,70,10 
10 0S»0CHECK(KU,0) 

IP ( PS-PCH2 ) 20,20,50 
20 TS*TSS( PS I 

CALL OENStl »TSfZEf ZEfDSLtOSVf 1 ) 
IF C DS-OSL ) 30,60,40 
30 IF ( OS— DSV I 50 , 60 1 60 
40 KR*2 

TS » TS+.Ol 
GO TO 80 
50 KR*3 
TS-1.2 
GO TO 80 


REGION 1 

60 KR»1 

GO TO 110 
70 T S*TSS I PS I 
GO TO 110 

REGIONS 2 ANO 3 

80 CALL QMUST(DS) 

CALL OOTIDS, TS) 

CALL Q2DTfDS,TS> 
TS*SOLVE(TS tTSFfOTSF ) 

VERIFY REGION 


IF I PS-PCH2) 1 10 t 110*90 
90 IF ITS-TCH2) 110,100,100 
100 KR»2 

110 T«TS*TCONV(KU> 

RETURN 

ENO 
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noonnoon 


«IBFTC TSFI 

FUNCTION TSFITSl 
C 

c 
c 
c 
c 


VERSION MARCH 


FUNCTION USEO TO SOIVE FOR TEMPERATURE TS GIVEN PRESSURE 
ANO DENSITY 

COMMON /PRHOTT / PS.O .T 

COMMON /CCNSTS/ TAUC.RHOA.RHCB.TAUA.E *R 

P ST AT F«100oi*R-0 /TS*11.*0 *IUCALCCTS> *0 *ODTAITS»)» 

TSF-PSTATE-PS 

RETURN 

ENTRY OTSFI TS > 

01 SF«R*0*l I [!o>0*0«UOT* its . .O.OCALCIIS ll-IS WI0.020IAIIS . 
l+QTDCTS 

OTSF*D1SM(-1000./iIS *T S )) 

TSf«OTSF 

RETURN 

FRO 
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S1BFTC ENTH1 

SUBROUTINE ENTHI KU f TT ,O f H) 


-VERSION MARCH 1,1972- 


THIS ROUTINE COMPUTES ENTHALPY GIVEN THE TEMPERATURE PARAMETER TT 
Iwn THE DENSITY D. I/O UNITS ARE SPECIFIED BY KU. 

IF SATURATION VALUES ARE NEEDED. THIS ROUTINE MUST BE CALLED 
WITH DL AND OV INPUT AS D. 

ENTHALPY IS RETURNED IN H. 


COMMON/ IERROR/ IROUT 

C0MM0N/CCNV6/HC0NVI 5 ) 

COMMON/SICOF / PSI1. PSI2.PSI3.PSIA. PSI5 
COMMON /CONSTS/ TAUC.RHOA.RHOB.TAUA.E ,R 
I ROUT-4 

TS»TCHECK«KU,KR.TT> 

OS-OCHECK(KU.O) 

CALL QMUST IDS) 

CALL QMUST 2 ITS I 

T PSI0-*(PSI34T+PSI2I*T4PSI1*(PSIA4PSI5*TI*AL0G(T) 

”sll»MIVT .PSI5.I1..AL0GII I. 

!lUl500^lmiu.AI)S*<<ICAlCIIS|.TS.OTOItS.ADS.<.0 T l0S. T Sm 

H* (H1+H2) *HCGNV(KU ) 


RETURN 

ENO 
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SI BFTC ENTl 


SUBROUTINE ENTtICU, TT, 0*$ > 

VERSION MARCH It 1972 

THIS ROUTINE COMPUTES ENTROPY GIVEN THE TEMPERATURE PARAMETER TT 
AND THE DENSITY 0. I/O UNITS ARE SPECIFIED BY KU* 

IF SATURATION VALUES ARE NEEOED* THI $ ROUTINE MUST BE CALLED TWICE 
WITH DL AND OV INPUT AS D« 

ENTROPY IS RETURNED IN $• 

CGNMON/SIICGF / PSUtP$I2*PSI3tPSI4, PSI5 
COMMON /CONSTS/ TAUC « RHOA, RHOB» TAUA t E ,R 
COMMON/ I ERROR/ I ROUT 
COMMON/C€NV4/$CONV(5 I 
I R0UT*5 

T S*T CHECK ( KU«KR« TT ) 

DS*DCHECK<KU»0> 

CALL QMUST(DS) 

CALL 0MUST2CTS I 
T-IOOO./TS 

PSIT«2**P$I3*T *PSI2*PSI4/T ♦PSIS*! 1*+AL0G< T i) 

SSS»-R*f ALOGI DS I>DS^ C QCALC ( TS )-TS*QTD(TS! I l-PSIT 
$»SSS*SCCNV( KU ) 

RETURN 

END 
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SIBFTC TEMPP1 


SUBROUTINE TEMPPHIKU, P,H,T,O f OL,OV,KR ) 
c VERSION MARCH 1,1972 

COMMON /CONV1/OCONVI 5 ) 

COMMON /C0NV2/TCGNV( 5 ) 

COMMCN/CC NV6/HCONV ( 5 ) 

COM MON/PHC ALL/PS , HS ,SS , 

COMMON/CHECKS/OCHl,DCH2,PCHl,PCH2,PCH3,TCHl,TCH2,TCH3,DST,TST » 

lHSCHlf HSCH2 
COMMON/ IERROR/ I ROUT 
EXTERNAL TSHF 
PS*PCHECK <KU,KR,P) 

IROUT-6 

HS*F/HCGNV ( KU ) 

IF (HS-HSCHU 20,10,10 
10 IF (HS-HSCH2) AO, 40, 20 
C 

C INPUT H - OUT OF RANGE TAG 

20 WRITE(6,30U HS , HSCH1, HSCH2 

301 FORMAT (10H0INPUT H * ,G14.6, 29HJ/G IS OUT OF RANGE OF HMIN= 

1 , FA* 1 , 10HAND HMAX * ,F7.1, 3HJ/G ) 

C 

AO IF (PS-PCH2* 140,140,130 
130 TS1*TCH1 
TS2-TCH3 
GO TO 110 
140 TS-0.0 

CALL DENS(1,TS,PS,ZE,DL,CV,1> 

IF CKR-1) 50 » 70 1 30 
50 CALL ENTH(1,TS,DL,HSL ) 

CALL ENTH(1,TS,0V,HSVI 
IF (HS-HSU 90,70,60 
60 IF (HS-HSV) 70,70,100 
C 

C REGION 1 

C 

70 KR-1 

80 CALL 0ENS(1,TS,ZE,ZE,CSL,DSV,1) 

OL*DSL*OCCNV ( KU ) 

DV*OSV*OCON V ( KU ) 

GO TO 120 
C 

C REGION 2 

C 

90 KR*2 

TS1*TCH3 
PS*PS*1.000 1 l 

T S2* 1000* /( 1000. /TS-l.E-5) 

GO TO no 
C 

C REGION 3 

C 

100 KR*3 

TSWOOO./QGOO./TS+1.E-5) 

PS* PS*. 99988 
TS2*TCH1 
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26 
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28 

29 
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REGIONS 2 AND 3 

110 CALL ROOTX (TSl*TS2*HSfTSHF*TS) 
CALL DEN$(l,TStP$,DS*ZE*ZE,KR) 
0*DSADCCNV(KU) 

VERIFY REGION 

IF ( PS-PCH2 I 120*120 *150 
150 IF C TS-TCH2 ) 170*170*160 
160 KR«2 

GO TO 120 
170 KR*3 

120 T»TS*TCCNV(KU) 

RETURN 

END 
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$ IBFTC STEMPS 


SUBROUTINE TEHPPS (KUt Pt S, T # 0,DL,DV,KR ) 

VERSION MARCH 1, 1972— 

COMMON /CONVl/CCONVi 5 ) 

CGMM0N/CCNV4/ SC0NVI5) 

COMMON /CCNV2/TC0NV( 5) 

CCMMON/PHC ALL/PS * HS *SS 

COM MON/CHECKS/ DCH1 * DCH2# PCH1 » PCH2,PCH3 f TCHl f TCH2 f TCH3 t DST*TST f 
1HSCH1 *HSCH2 
CCMMON/IERROR/I ROUT 
EXTERNAL TPSF 
IR0UT*7 
SMAX*1 3* 26 
PS*P CHECK ( KU # KRt P ) 

SS = S/SCCNV ( KU ) 

IF (SS .LT • O.O) GO TC 20 
IF ( SS.IE.SMAXI GO TO 40 


C 

C INPUT S - OUT OF RANGE TAG 

20 VIR I TE( 6f 30 1 ) SS,SMAX 

301 FORMAT (iOHOINPUT S * f G14.6t 47HJ/G-K IS OUT OF RANGE OF 
10 AND SMAX* tF7.lt 5FJ/G-K ) 

40 IF (PS-PCH2) 140 1 140 1 130 
130 TSlxTCHl 
TS2*TCH3 
GO TO 110 
140 TS*0 .0 

CALL DENS(l,TS,PStZE f CL f DV,l) 

IF (KR-1) 50 » 70 # 50 
50 CALL ENT( ltTS,OLtSSL ) 

CALL ENT ( 1 1 TS, OV » SSV ) 

IF ( SS-SSL ) 90 1 70 » 60 

60 IF ( SS-SSV) 70 f 70» 100 


SMI N*0 • 


REGION 1 


70 KR* 1 

80 CALL OENS(ltTS,ZEtZE,CSLtDSVtl) 
DL*DSL*DCONV ( KU ) 

OV*DSV*DCGNV ( KU ) 

GO TO 120 

REGION 2 


90 KR*2 

TSl*TCH3 
PS*PS*1. 00011 
TS2*TS*1. 00001 
GO TO 110 

REGION 3 

100 KR*3 

TS1*TS*. 99999 
PS* PS*. 59988 
T S2* TCH1 
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REGIONS 2 AND 3 

110 CALL ROOT X ( TS1 » TS2, SS, TPSF , TS ) 
CALL DENS(1, TS,PS,DS,ZE, ZE,KR) 
0*DS ♦DCCNV ( KU ) 

VERIFY REGION 

IF (PS-PCH2) 120,120,150 
150 IF ( TS-TCH2I 170,170,160 
160 KR*2 

GO TO 120 
170 KR*3 

120 T»TS*TCCNV(KU) 

RETURN 

ENO 
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SI6ETC TSFfl 

FUNCTION T SHF I TS 1 

C VERSION MARCH 1,1972 

CON PON /PNC ALL / PS ,HS ,S S 
KR*0 

CALL 06NSI 1,TS,PS,0S,/E,2E,KRI 
CALL ENTM 1*TS*CS#HSC> 

T SHF*HSC 
RETURN 
f 

ENTRY TPSFI TS I 
KR * 0 

CALL DENS ( 1, TS, PS,0$ ,2E ,2E , KR 1 
CALL FNTI 1*TS,DS«SSC> 

TPSF * S SC 
T$HF*TPSF 
RETURN 
FNC 
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$ IBFTC CPPRL1 

SUBRCUT INF CPPRL IKU*T«OtCP tC V* GAMMA* C) 

t VERSION MARCH it 1972 

C f 

C THIS SUBROUTINE RETURNS THE FGLLOWI KG TO WASP IK USERS UNITS* 

C SPECIFIC HEAT AT CONSTANT PRESSURE *CP 

C SPECIFIC HEAT AT CONSTANT VOLUME «CV 

C SPECIFIC HEAT RATIO *GAMMA 

C SONIC VELOCITY * C 

C 

C THE PART I AL S PTV AND POT EXPLAINED BELOW ARE RETURNEO IN COMMON* 

l 

COMMON /COC PO/ COC 1 «COC2*COC3 . 

COMMON /S I COF/ Cl tC2*C3tC4,C5 

COMMON /PARTLS/ PTVtPOT 

COMMON /CON STS/ TAUC «RHOA tRHOB * TAUA *E «R 

COMMON /CON V4/SC0N VI SI 

COMMCN/COKV 5/CCOK Vi 5 1 

COMMON/ IERROR/ I ROUT 

IRCUT*8 

T 5*TCHECK ( KU*KR«T I 
TT* 1000*/T S 
OS* DCHECK I XU* 0 I 
CALL OMUSTIOSI 
CALL OMUST 2ITS I 
CQ2T20*Q2T2C< TS I 
COC ALC*OCALC4 TS I 
COJOT*QOT I OS tTS I 
CCJOaTDdS) 

CQ2OT*Q2CTI0S*TS> 

CQ202T*02C2TI0StTSI 

CV * -2 **C3*TT +C4/TT-C 5— R*DS*TS*TS*C 02T20 

C PTV IS PARTIAL OF P BY T INOT TAUI 

C~- POT IS PARTIAL OF P BY RHO 

PTV«R*OS*U«+OS*< COCALC+OS*COOT-TSA(CQTD+DS*C02DTIII 
POT»R*TT* ( 1 **DS* ( 2 **CQC ALC +OS* 1 4**C0DT +0S*CQ202T 1 1 I 
OHOT*-i.*C3*TT+C4/TT-C5*R*< 1 . ♦DS’M COCALC+D S*CQOT-T S* ( OS*CQ2DT+ 
1CGTC+T S*C02T20 I 1 1 

DHOO~R* t TTACQC ALC *1000 **COTO*D$*(TT*( 3. *C00T+DS*C0202T 1+1000. * 
1C020T) I 

CP * DFOT-DHCCMPTV/POTI 

GAMMA*CP/C V 
CP*CP ASCONV (KU I 
CV*CV«SCONV4KUI 
GAMMAP*GAMMA * 10** POT 
CS«0*O 

IF ( GAMMA P*GT • 0*01 CS* 1000 .♦ SORT IGAMMAPI 
C« CS*CCCKV ( KU ! 

RETURN 

END 
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fIBFTC SRTNSB 


SUBROUTINE SURF ( KU*KR * T IN* SURFT ) 

* VERSION MARCH i, 1972 

THIS ROUTINE CALCULATES THE SURFACE TENSION OF LIQUIO MATER AND 
THE LAPLACE CONSTANT 

COP MON/ I ERROR/ I ROUT 
COPPON/CCNV9/STCONVI 5 ) 

COMMON /LAPLAC/ ALC 
DIMENSION A( 5 ) 

DIMENSION X(5) 

DATA IAm*I»lt5)tB f TK /0. 11609368 , 1*1214067 E-3» -5.7528052 

IE-6, 1.2862744 E-8, -1.1497193 E-ll* 0.83, 647.30 / 

C UNITS OF G - N/S**2 

DATA 0 / 9.80665 / 

I ROUT-11 

C T IS OEG K 

TAU-TCHECK C KU.KR , T IN 1 
T - 1000. /TAU 
SURFT-0.0 
ALC-0.0 

IF (T.GT.TKI RETURN 
X(l>» TK-T 
XI2I- XllltXIll 
XI3I- X(2)*X<1) 

XI4I- X(3)*X<1) 

XI 5 )- X(4)*xm 
Y « ( All ) *X(2) )/( 1.0+B*X (1 ) ) 

00 1 N-2,5 

1 Y - Y+AINI * X(N» 

SURFT - Y 

C UNITS OF SURFT MUST BE DYNE/CM 

C UNITS OF ALC IS MM. 

TR- T/647.30 

IF (TR.GT. .998) GO TO 2 

CALL DENS (KU,TAU, ZE, ZE* 0L*DV«1 ) 

ALC • SORT (SURFT/ ( G* ABS (DL-DV)MOOO. I I 
C— CONVERSION FACTOR FOR RESULTS TO BE IN MM AS IN THE TABLES 

2 ALC » ALC * 31.622777 
SURFT-SURFT*STCONV(KU ) 

RETURN 

END 
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SIBFTC OTH 


FUNCTION DTHERMI XLM) 

FUNCTION USED TO SOLVE EQ«(B52l FOR THERMAL CONDUCTIVITY 


COMMON/ ITER AT/TRf TO* Tl*T 2 fT 3 *T 4 ,T 5 *T 6 *T 7 *T 8 
DOUBLE PRECISION TO* T 1 *T 2 * T 3 *T 4 * T 5 *T 6 *T 7 , T 8 

TCALC* ( XLM* ( XLM*( XLM* ( XLM*( XLM* ( XLM*( XLM*T 8 +T 7 )*T 6 )+T 5 )*T 4 )*T 3 I* 


IT 2 )+T 1 )*XLM*T 0 
OTHERM*TC ALC-TR 
RETURN 

ENTRY DOTH ( XLM ) 


DERIVATIVE USED TO SOLVE EQ.(B52> IN NEWTON RAPHSON ITERATION 


C 0 TH«<XLM*tXLM*tXLM*<XLM*(XlM*{XLM* 8 .*T 8 V 7 .*T 7 >* 6 .*T 6 )* 5 .*T 5 >+ 

14**T4)*3«*T3)+2«*T2)*XLM+T1 

DTHERM*DDTH 

RETURN 

ENO 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
19 


66 



ooonooono 


SIBFTC STHCND 


SUBROUTINE THERM (XU* KR f P IN, T IN, DIN f EXCESK, TCOND ) 

VERSION MARCH 1,1972 

SUBROUTINE CALCULATES THE THERMAL CONDUCTIVITY IN INTERNAL 
UNITS OF W/CM-K AND CONVERTS TO USERS UNITS 
EQUATIONS ARE THE INTERNATIONALLY AGREEO UPON ONES IN REGIONS 
NHERE SAME ARE AVAILABLE AND ARE PROPOSED EQUATIONS IN OTHER 
REGICNS. 

THE NEAR SUBCRITICAL REGION IS THE AUTHORS FIT 

t 

C0MM0N/CCNV8/KC0NV<5 ) 

COMMON/ IERROR/ I ROUT 

COMMON/ ITERAT/TR, TO, T1,T2,T3,T4,T5,T6,T7,T8 
REAL KCCNV 

COMMON /CRIT/ RHOCRT , PCRT, TCRT 
CCMMGN/TPARAM/TAU 
EXTERNAL DTHERM * DOTH 
DIMENSION CFC( 5* 2 ) 

OAT A CFC/- .5786 1540* 1 .45 746404, * 17006978* • 13348045* • 3278399 IE-1 * 
1-. 70 859254 *.94131399, .64264434E-01 , 1.85363188,1.98065901/ 

DOUBLE PRECISION PROP 

OOUBLE PRECISION A5 C 5 ) * B5 ( 4 ) ,C5( 4 ) * T1 * SUH1 , SUM2 , SUM3 
OATA AO, A1,A2, A3, A4/-. 17384732,. 82350372, -1.55213983, -.12626 138, 

1 2.83922425 / 

OATA A5 / -.9224700000 * 6.728934102 , -10.11230521, 6.996953832* 

1 -2. 316062510 / , B5 / -.2095427600, 1.320227345, 

2 -2.485904388, 1.517081933 / , 

3 C5 / .08 104183 1 A 7 , -.4513858027 , .8057261332, -.4668315566 / 

OOUBLE PRECISION A10 ( 2 I ,B 10 ( 2 ) ,C10< 3 ) , 010 ( 6 ) , A,B,C , T1 1 * T22 , T33 
OATA A10 / .01012472978, .05141900883 /,B10 /663742.6916, 

1 1.388806409 / , CIO / 338855.7874, 576.8000000, .2060000000 /, 

2 010 / .000002100200454, 23.94090099, 3.458000000,13.63235390 , 

3 .01360000000, .C07852600000 / 

DOUBLE PRECISION A8 ( 9 ) ,B8 ( 9 ) , C8 , E9( 3 > , TO, T2, T3, T4, T5, T6, T7, T8 
OATA A8 / 1.365350409, -4.802941449, 23.60292291, -51.44066584, 

1 38.86072609, 33.47617334,-101.0369288, 101.2258396, 

2 -45.69066893 / , 

3 B8 / 1.514476538, -19.58487269, 113.6782784, -327.0035653, 

4 397.3645617, 96.82365169,-703.0682926, 542.9942625, 

5 -85.66878481 / , C8 / 1.017179024 / 

DATA E9 / 50.60225796 , -105.6677634 , 55.96905687 / 

I RQUT*10 

PMN*PCHECK ( KU,KR, PIN ) 

TAU*TCHECK (KU,KR,T IN ) 

DS*DCHECK (KU,OIN ) 

C CONVERT TAU AND PMN TO VARIOUS UNITS 

TK-1000./TAU 
TR * TK/TCRT 
T * TK-273.15 
PBAR * PMN*10. 

PR*PMN/PCRT 

C CUT OF RANGE CHECK ON PRESS ANC TEMP. 

IF ( PBAR. LT • 1.0 .OR. PBAR. GT. 500 . \ WRITE(6,151) TIN, PIN 
IF (T.LT.O.Q.OR.T.GT.7GO.) WRITE(6,151) TIN, PIN 
151 FORPAT ( 1 HO , 5H T *,F12.4,8H OR P =,F12.4, 64HIS OUT OF RANGE, 
1RETURNE0 THERMAL CONDUCTIVITY IS EXTRAPOLATEO ) 

C 

C CHECK FOR REGION I 

C 

IF (T. LE. 350.. AND. DS.GT. RHOCRT) GO TO 100 
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CHECK FOR JAGGED LOWER BOUNDARY OF REGION IIIOR UPPF.F. PART 62 

OF REGION II 63 

66 

IF < PBAR. GT. 650.. AND. T.LT. 550.1 GO TO 80 65 

IF (PBAR.GT.35O..AN0.T.LT.500.) GO TO 80 66 

IF (PBAR. GT. 275.. AND. T.LT. 650.1 GO TO 80 67 

IF ( PBAR.GT. 225. .AND.T .LT.625.) GO TO 80 68 

IF (PBAR.GT. 175. .AN0.T.LT.600.) GO TO 80 69 

70 

EQUATION ( B60 I FOR P-1.0 BARS 71 

72 

10 ¥1 * <17.6 ♦.0587*T*.000106*T*T -6.51E-08*T*T*T)/1000. 73 

IF (PBAR.GT. 1.0001 GO TO 20 76 

TC0ND-V1 *KC0NV ( KU I 75 

GO TO 500 76 

77 

EQUATION (B66) FOR REGION II. 78 

20 ANS* ( ( 103.51+.6198*T-2 »771E-05*T*T I*DS+2.16821E*16/(T**6.2 )*DS*DS 79 

U/1000.+V1 80 

TCONO*ANS*KCONV(KU I 81 

00 TO 500 82 

83 

REGION I CALCULATIONS 86 

85 

100 CALL PSSS(PS) 86 

PREOD « (PMN-PSI/PCRT 87 

SUNl-I < < A5 < 5 I *TR+A5< 6 I ) *TR*A5 < 3 i ) *TR+A5< 2 1 )*TR*A5 (1) 88 

SUM2»< <B5(6)*TR+B5<3n*TR+B5<2) )*TR+B5<1I 89 

SUM3»< (C5 ( 6 )*TR+C5( 3 1 I *TR+C5( 2 > )*TR*C5< 1 1 90 

TCONO* ( SUNl ♦ (SUH3*PREDD + SUM2)*PREDD )*KCONV(KU) 91 

GO TO 500 92 

93 

CHECK FOR REGION IlI-USING BOUNDARY EQUATION (B56) WHICH 01VIDES 96 

REGIONS III AND IV. 95 

THEN SEPARATE HATCHEO REGION WHERE NO EQUATION EXISTS FROM REMAINDER 96 

OF REGION IV 97 

80 I EQUA-10 98 

I F( T.GT .650. ) GO TO 300 99 

PB0UND«E9< 1)+E9(2)*TR+E9(3)*TR*TR 100 

IF ( PR.LT.PBOUND I GO TO 300 101 

I F ( TK. LT. TCRT. ANC.OS • LT.RHOCRT ) GO TO 600 102 

IEQUA-8 103 

106 

EQUATION <B52> IS SOLVED BY ITERATION 105 

106 

200 NTR-0 107 

PROP-PR 108 

C0N-PRDP-C8 109 

TO - A8 ( 1 )♦ C0N*B8 ( 1 ) 110 

T1 » A8 ( 2 I ♦ C0N*B8 ( 2 ) 111 

T2 * A8 ( 3 )♦ C0N*B8 < 3 ) 112 

T3 - A8( 6 C0N6B8 ( 6 ) 113 

T6 - A8 ( 5) ♦ C0N*B8<5) 116 

T5 = A8(6)+ C0N*B8 ( 6 ) 115 

T6 - A8(7t« C0N*B8<7) 116 

T7 * Afi ( 8 )♦ C0N*B8 ( 8 ) 117 

T8 - A8<9|* C0N*B8 ( 9 ) 118 

C OSE CONDUCTIVITY BASEC ON BOUNOARY AS INITIAL ESTIMATE 119 

PPR - PR 120 

PR-PBOUNO 121 

GO TO 300 122 
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210 PR 3 PPR 
XHI 3 .55 

X 3 ( ANS*3.+XHI*2. >/5. 

IF( PR.LE. 1.05. AND. TR.LE. 1.05 ) X*ANS+.005 
TCOND 3 SOLVE (X»DTHERM*CDTH) *KCONV ( KU ) 

THIS EQUATION OOES NOT ALWAYS CONVERGE NEAR THE BOUNDARY WHERE 
IT SHOULD. SWITCH TO AUTHORS EXTRAPOLATION IF THIS HAPPENS. 

IF (TCCND.LE.O.OI GO TO 400 
GO TO 500 
300 DON 3 1.00+0 

B-CB10<1) *PR**1«63) / I C0N+B10 ( 2 ) *PR**3.26) 

C* (ClO(l) *PR*#1 .5 ♦ C10I2M / B - C10I3I 
CSP 3 C 

TEST * DCN-B*D10(1)/TR**7 
T11=(A10( 1 )*PR+A10(2> )*TR**1.445/TEST**CSP 
T22= D10(2)*PR**4 *EXP(-9.0*D10( 3 )*( TR-1.0) ) / ( D0N+D10 ( 41 
1 /PR**12 ) 

T33* 010(5) - 010(6) *PR *EXP(-D10(3) * (TR-1.0)) 

ANS 3 Til + T22*T33 
IF (IEQUA.EQ.8) GO TO 210 
TCOND*ANS*KCONV(KU) 

GO TO 500 

HATCHEO REGION WHERE NO EQUATION EXISTS IN THE REFERENCES. 

C AUTHORS OWN EQ USED HERE WITH 1 BAR EQ. 

400 VI « (17.6 +.0587*T+.000104*T*T -4.51E-08*T*T*T)/1000. 

XX- AL0G10 ( DS/RHOCRT ) 

KJ«1 

IF (XX. GT.-. 39794) KJ»2 

**CFC ( 1 1 X J )+ ( ( (CFC(5»KJ)*XX+CFC(4tKJ) ) 4XX+CFC( 3, K J ) )*XX 
1 +CFC(2,X J) )*XX 
TC0ND 3 (10.**Y+V1)4KC0NV(KU) 

500 CONTINUE 
C 

C REACTING CONDUCTIVITY IN THE NEAR CRITICAL REGION BY SENGERS 

DRHOC 3 ABS ( OS - RHOCRT ) / RHOCRT 
DELAMB=0. 

I F ( DRHOC. GT. .6) GO TO 520 
DELTC 3 ABS (TR-1.) 

RAT = DS/RHCCRT 

IF (ORHOC.LT. l.E-4) GC TO 510 
IF (DELTC. LT.l.E-7) GC TO 502 
XBETA 3 DELTC**. 35/DRFOC 
IF (XBETA. GT.. 4) GO TO 506 
502 DELAMB* 11.6E-5 / ( SQRT < RAT ) *0RH0C**1.71 ) 

GC TO 520 

506 IF ( XBETA. GT .3. ) GO TC 510 
XB* ALOGIO ( XBETA ) 

ARAT 3 ( < ( A4*XB+A3)*XB+A2)*XB+A1)*XB+A0 
DELAMB* 11 .6E-5/ ( SQRT ( RAT )*DELTC**.6)*10.**ARAT 
GC TO 520 

510 IF (OELTC. LT.l.E-7) GO TO 525 

0ELAM8 = 11.6E-5/(SQRT(RAT)*DELTC**.6) 

520 EXCESK 3 DELAMB 
RETURN 

525 EXCESK 3 l.E30 
RETURN 
END 
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SI8FTC C Y V f SC 

SUEHCUTINf VISCIKU*KR.TIN,PIN«CIN« S V l SC ) 


r — VERSION MARCH 1,1*72 

c 

C CALCULATE THE VISCOSITY GIVEN TAUfP.ANU 0 IN USER-S UNITS MJ. 

C ANSWER RETURNED IN USER-UNITS IN SVISC. 

f 


CCMPON/ IERRCR / IROLT 
COMMON /CR IT/ RHCCRT ,PCRT , TCRT 
COMMON /T PAR AM / TIN 
COMMON /CCNV7/MCGN VI 5* 

REAL MCCNV 

D I MENS It N A I 5 ) , £ ( 3 ) . C I 31,0(31 ,COE (5,2) 

OATA (AID, 1*1 • 51 /241 .4000,0.38282095,0.21628302, 0*1498694, 

1 0.47118801 i , 

2 mi>'I-l,3> /263. 4511, 0*4219836, 80.4000 / , 

3 (C( I ) , I * 1 , 3) / 586.11587, 1204.75394, 0*4219836 /, 

A (D(l>, 1 = 1,3) / 111.35647, 67. 32080 1, 3.2051470 / 

CAT A C0F/-6 .4556581 . 1.39494 36, .302 590 83, .10960682, • 1523003 1E-0 1 , 
1-6 .4608381 , 1.616332 1 0 , 1 . C 7C 9 7 70 5 , - 1 3 . 93 8 , 30 . 1 19832 / 

18001*9 

C 

c TK IS DEG K, T IS DEG C, TR IS REDUCED TEMP 

TK*1000./TCHECK(KU,KR,TIN ) 

T * TK- 273.15 
TR*TK/TCRT 

C P IS BARS, PMN IS MEGA NEMTCNS/M*M ,PR IS REDUCED PRESSURE 

PMN*P CHECK ( KU , K R ,PIN! 

P =* 10.04PMN 
PR*F MN/PCRT 

C DC IS G/CC, SPVR IS REDUCED SPECIFIC VCLUHE. 

OD=OCHECK I KU, D I N> 

SPVR* RHOCRT/DO 
C 

C CHECK FOR TUT OE RANGE ON P ANO T 
C 

IF I P .LT. .99 .OR. P .GT. 800.01) URITE16,i01) T.P 
IF I T .GT. 800.0 .OR. T .LT. 0*0 ) MRITEC6.101) T,P 

101 FORMAT I 1H , 48H OUT OF RANGE. AtiSkER IS EXTRAPOLATED FOR T« 
1.F12.4* 4H P* , F 1 2*4 ) 

IF (CO.GT .RHOCRT. AND. T *L T .300. ) GO TO 100 
GO TO 110 
C 

C REGION I 

C FOR TEMP — 0 TO 300 CENT ANO PRES — PSAT TO 800 BAR 

C 

IOC XI * 10.0 **IA(2)/(TR-A(3) ) ) 

CALL PSSSIPSS) 

PSR*PSS/PCRT 

X? = 1.0 ♦ ( PR— PSR)*AI4)* (TR— A( 5) ) 

SVISC * (All) ♦ X 1 *X2 )/10.0**6 ♦MCONV I KU) 

RETURN 

C 

C CALCULATE VISCOSITY FOR 1.0 BAR NEEDED FOR REGIONS II,III«IV 

C 

110 VISC1=(P( ) >*(Tfi-B(2) >♦8(3)) 

IF (P.LE.I.O) GO TO 1000 

c 

C CFFCK FCfi TEMPERATURE RANGE 8HERE NO CURVE EXISTS 

C 

IF IT .GE.300..AND.T.LT .374.15) GO TO 400 
IF (Oi).LT • RHOCRT .AND*T .LT.300. ) GO TO 200 
GC TO 300 
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C 

c REGION II 

C— FOR PRES — 1 TO PSAT BAR AND 100 TO 300 CENT 
C 

200 SVISC « ( V I SC 1 - 1.0/$PVR*(C( 1 )-C 1 2 1* I TR— C 131) II /l 0«**6*MC0NV C KU1 
RETURN 


REGICN MI 

-~FCR PRES — 1 TO 800 BAR AND 375 TO 800 CENT 

300 SVISC -IVISC1 ♦ DI11/SPVR +DI 2 1 / I SPVR* SPVR 1 ♦ 0( 31 / ( SPVR*SPVR*SPVR 
1 1 1/10. 0**6 **4CGNV(KU> 

RETURN 

C 

C AUTHORS EXTRAPOLAT ICN USEC FOR REGICN IV 
C 

400 X-ALCG10I 1 ./SPVR) 

KJ«1 

IF (X .GT.-. 124938731 KJ*2 

V* X*< X*(X*IX*CCF(5«KJ l+COF (4«KJ 11 +COF (3 *K J! 1*C0F < 2,KJ 1 1 *COF ( 1 .KJ > 
SV I SC* I V I SCI/ 1 .E6+10 •**{ Y*l. I /. 01921 *MCONV (KU> 

RETURN 

100C SVf SC-VISC1/10.**6*MCGNV(KU1 
RETURN 
ENC 
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APPENDIX F 


TEST PROGRAM WITH OUTPUT 

The following tables have been generated by WASP to facilitate comparing results to 
the ASME Tables (ref. 1) and the International Skeleton Tables (refs. 1 and 2). No at- 
tempt was made to reproduce the entire reference tables; only a select number of points 
were chosen at even intervals representative of the total range. The values in the fol- 
lowing tables are in the same units and similar form as the reference tables. Results of 
comparisons between the calculated and the tabulated values are discussed in the main 
part of the text. 


$ IBFTC MTNASP 


C 

C 

C 

C 

C 

C 

C 

C 


cc»«»*on/phoi»ty/ku,ol,dv,hl,mv,s.sl.sv.cv,cvl.cvv,cp,cpl,cpv.gamma, 

lGANMALtGANNAV.C,CL.CVP,NU,NUL,NUV,K,KL.KV,SIGMA,EXCL.EXCV,EXCESK 

DIMENSION PSIAI 12 I* TF 1 11 )• VOL< 12» 11) tHOUTIl2.il> .SOUTI 12* 11 1 
DIMENSION P6ARTCI13) » TCOUT C 20» 13) tTC0UT2(3l» 13 ) 

DIMENSION T( 200) »PI200 ) « VLI 200 ) t HV0UT 1 200) »SV0UTI 200) »HL0UT I 200 ) 

1 1 SLOUT ( 200 ) » VV ( 200 ) 

DIMENSION PNEARI 12)»TNEA«(6) 

COMMON/LA PL AC/ ALC 

115) 


MASTER TEST PROGRAM FOR WATER AND STEAM PROPERTY PACKAGE 


PART 1 


COMPARE SATURATION PROPERTIES 
AS A FUNCTION OF TEMPERATURE 


OF ASME STEAM TABLES PAGES 83-88 
35-705F IN INCREMENTS OF 10F 


KU*3 
NPT* 1 
Tl*494. 

10 KR*1 
P1*0.0 

CALL WASP(lt3,Tl t Pl,RH0,H,KR) 
T(NPT)*T 1-460* 

P(NPT) *P1 
VL(NPT)*1./0L 
VV (NPT 1*1 • /DV 
SLOUT ( NPT I *SL 
SVOUT ( NPT )*SV 
HLOUT ( NPT ) *HL 
HVOUT ( NPT ) *HV 
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C 

C 


NPT-NPT+1 

T1*T1+12.0 

IF (T1.LE.1165. ) 60 TC 

npt*npt-i 


10 


VL 

SV 


FT3/L8H 


VV 

) 


HL 


BTU/L 


C 

c 

c 

c 

c 

c 


PRINT SATURATION RESULTS 
NR I TE ( 6 V 1 ) 

1 VRITE(6l2| ,5lH COMPARA0LE T0 ASME TAB LE NO. 1 PAGES 83-88 ) 

M I N*0 

00 20 J*1,NPT 
J J*NPT-J«-1 
NLIN-NLIN+1 

USVOUTUJ) T,JJ,,PUJ, » VU JJ >» vv *JJ*»HLOUTIJJ> f HVOUTUJ),SLOUT(JJ) 

IF INLIN.LT.50) GO TC 20 

NLIN-0 

MRI TE(6,1 ) 

MR I TE (6,2 ) 

2 FORMAT ( 1H0,90H T-F P-PSIA 

1BB HV SL BTU/LBH-R 

20 CCNTInIe” F5 - 0 * F10 *^* 2f: »2.6.2F12.3,2F12.-%* 

WRITE(6,21I 

21 FORMAT C 1HI ) 

PART 2 ASME TABLE NO. 3 PAGES 97-203 

PROPERTIES OF SUPERHEATED STEAM AND COMPRESSED MATER 
TABLE IS 2 PAGES 32-750F AND 750-1500F 
FOR EACH SET OF ISOBARS 

MILL COMPARE 12 ISOBARS FOR VARIOUS TEMPERATURES 

Kai5So!5 /1 * 0,5 * 0,25 * ,loc * ,200 * ,500 *’ looo ** l50<) *» 2000 ** 5000 -' 100 o 0 

DO T 50 T I-iai 50 * ,100 * ,150 * ,300 * ,500 * ,700 * ,900 * ,ll00 * ,1300 * ,l500 * / 

TIN*TFUI*460. 

CO AS J*I,12 

KR»0 

KU»3 

CALL MASP(l,3,TIN,PSIA(J),D f H,KR) 

IF (KR.EQ.l) GO TO 45 
VOL (Jtl ) * 1 • /O 
HOUT(J,II*H 
SOUT ( J , I ) *S 
45 CONTINUE 
50 CONTINUE 

MRITE(6,41 I 
MRITEI6.42I 
CO 60 J*l,12,3 

MR I TEC 6,44) P$I A ( J) , PS I A C J*1 J , PSI A( J42 1 
00 60 1*1,11 

WRITE! 6, 43) TF(I),VOL(J, I ) , HOUT I J, I ) , SOUTI J , 1 1 , VOL f J«l. I 1 HflllTi 1*1 

l.I».S0UTfj4ia»,*0UJ*2,Il f H0UTIji2 f n f SOUTCjI 2 .il J * l ’ n ’ H0UT,J+1 

41 FORMAT ( IH1 ,20X, 36H COMPARISON POINTS FOE TABLE NO. 3 1 

43 FO C S;A% T !!H°irr^^°, LUP f:^ THALPY » ENTROPY * F0R PRf SSURE LISTEO ) 

FORMAT ( 1 H »F6,0,3(F12.5 t FlU2 t Fi0.A,5Xn 
A 4 FOR PAT ( 1H0 ,6X*3( F10*0 f AHPSI A # 20X ) ) 

60 CONTINUE 

PART 3 ASME TABLE NC« A PAGES 208-220 

THF P rBlJ?rA? F ocr?n? HEATE0 STEAM AND COMPRESSED MATER IN 
tut UKiTICAL REGION 

n«c\'MQ"")“i:r 0 '' 3C ‘ 0 -' MOO - ,5l ‘ 0 -’” 00 -' 3260 -' 3 “ 0 -’ 3 '' 00 -' 


35 

36 

37 

38 

39 
AO 
A1 
A2 
A3 
AA 
A5 
A6 
A7 
A8 
A9 

50 

51 

52 

53 
5A 

55 

56 

57 

58 

59 

60 
61 
62 
63 

6 A 

65 

66 

67 

68 

69 

70 

71 

72 

73 

7 A 

75 

76 

77 

78 

79 

80 
81 
82 
83 

8 A 

85 

86 

87 

88 

89 

90 

91 

92 

93 
9A 

95 

96 

97 

98 

99 
100 
101 
102 
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o o n 


C 

C 


DATA T NEAR/ 6 50. ,680. » 7 10 . ♦ 74C . » 770 . . 800 . / 

CC 80 1*1,6 
T I N = TNE AR ( I ) +460 • 

CC 70 J=1 , 12 
KR = 0 

CALL WASP ( 1 13, TIN,PNEAR( J) ,D,H, KR * 

VOL ( J, 1 ) * 1 • / D 
IE (KR.EQ. 1 > GO TO 70 
HCUT ( J, I )*H 
SCOT ( J » I ) =S 
70 CONTINUE 

80 CC N T INUE 

81 FORMAT ( 1H 1 ,20X , 40H COMPARISON POINTS FOR ASME TABLE NO. 4 1 

CC 90 J*1 , 12 , 3 

82 FCRM AT^lHO , 20X, 47H (VOLUME, ENTHALPY, ENT ROPY) FOR PRESSURE LlSTtD ) 
WRITE (6, 44) PNEAR(J) ,PN£AR(J+l),PNEAR(J+2) 

MUTEll’ot TNEARm,VCt(J.n.HOUTU,U.SOUT(J,I),VOUJ*l,n.HOUT( 
lJ+1,1) ,SOUT( J+l,I),VCL(J+2,l),H0UT( J + 2, l> , SCUT (J + 2, I > 

90 CON I INUE 

PART 4 ASME TABLE NO. 9 PAGES 278-279 
DATA PCP/1.,*., 10., 30., 60#, ICO. ,200., 400., 1000., 30CC., 6000., lOOUO 

^ OAT A TCP/32 . , 300 700 1 100 ., 1500 ./ 

109 F0RMAT(1H1,20X,42H COMPARISON POINTS FOR ASME TABLE NO. 9 

III FORM AT (IH/ ,4HPSl A, 12 ( 4X, F6.0 )/lH0, 4HTEMP/, 3H F /> 

WRITE(6,111> PCP 
DC 120 1 = 1,5 
T I N=TCP( I 1+460. 

CD 110 J= 1 » 1 2 
KR - C 

CALL WASP < 1 ,4, T I N, PC P ( J ) »D »H, KR I 

CPARYI J, I )=CP 

11C CONTINUE . 

HRITEI6, 112) TCP(I). (CPARYI J, I ) . J=l,12> 

120 CONTINUE 

112 FCRM AT < 1H0 ,F5.0, 12F10 .3 ) 

PART 5 VISCOSITY CHECKOUT 
DATA PBAR/1 . , 50. ,200 . , 35C. , 5C0. , 650. , 800./ 

OAT A TCENT / 10., 50.. IOC., 120., 140., 160., 180., 200., 220., 240., 260., 

. 300. . 320. , 340 . , 360. ,400. , 500. ,600. , 700. / 

DATA PSVS/1.,2.,5.,10.,20., 50., 100., 200., 500., 1000., 2000., 5000., 

^DAT A* TFVS/ 1500^ » 1450 • , 1400. , 1350., 1300. ,1250., 1200. ,1 150. ,1100. , 

1 1050.. 10C0., 950., 900., 850., 800., 750., 700., 650., 600., 550., 500., 

2450. . 400 ..350. .300. .250. .2 00. .150. .100. .50. .32./ 

WRITE(6«5C0) 

KC= 1 

CC 300 J= 1 , 7 
PI N = P8AR ( J ) / 10. 

CO 300 1=1,20 

TIN = TCENT ( I 1+273.15 

K R = C 

CALL WASP(l»8,TlNtPINtC»F,KR) 

IF (KR.EQ. 1) GC TO 3C0 
VCUT { I 9 J ) = MU* l • E5 
3CC CONTINUE 

WR I T E ( 6 ? 30 1 ) 

WRITF(6,302) P B AR 
CC 320 I = 1 » 20 

32C WRITE(6,303) TCENT ( I ) , ( VCUT ( I » J ) t J = 1 »7 ) 

3C1 FCRM AT ( 1H l t 42H VISCOSITY TABLE-INTERNATIONAL bCOK 


103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
1 22 
l?3 

124 

125 

126 

127 

128 

129 

130 

131 
13? 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 
169 
1 70 


76 



noon 


c 


c 

c 

c 


3C2 FCRMATUH0,7(8X,F6.0,4HBARS I) 
3C3 FOR HAT 1 1H0 ,F5.0,7 ( FI 2 .2, 6X ) ) 

HR I TE ( 6,500 ) 

KU*3 

CO 350 J-1,15 
PIN-PSVSI J) 

CO 350 1*1*31 
TIN*TFVS( 1)4460.0 
KR*0 


CALL WASP(1,8,TIN,PIN,0,H,KR) 
IF (KR.EQ.l) 60 TO 350 


THIS CONVERSION GETS FRON UINTS*3 OF PROGRAM TO UNITS OF TABL 
V0UT2I I, J)*MU*1. E6*l. 488 1639/4. 7880258 

350 CONTINUE 

HR ITE ( 6*351 ) 

NR I TEI 6,352 ) PSVS 
CO 355 1*1,31 

HR I TEI 6, 353 1 TFVS 1 1 ) , I VOUT 21 I * J I * J*l, 15) 

355 CONTINUE 

351 FORMAT 1 1H1 *45H VI SCOS ITY-ASME TABLE NO. 10 PAGE 280 1 

352 FORMAT(1HO,5H PS I A ,2X , 15F8 .0, /, 5H0 T-F ) 

353 FORMAT I IH , 2X * F6.0, l 5F8. 2 ) 


PART 6 THERMAL CONCUCTIVITY CHECKOUT 


360 CONTINUE 

OAT A PBARTC/1.* 10., 25., 50. *100., 150., 200., 250., 300. *350.* 400., 450. 
1 ,500. / 

WRI TE (6,500 ) 

KU* 1 

DO 400 J*1 ,1 3 
PIN=PBARTC( J)/10. 

DO 400 1*1,20 
TIN*TCENT( I)*273.15 
KR*0 


CALL WASP (1,16, TIN, PIN, 0,H,KR) 

IF (KR.EQ.U GO TO 400 
TC0UT(I,J)*K*l.E4 

400 CONTINUE 
WRI TE(6,401 ) 

WR I TE ( 6,402 ) PBARTC 
00 420 1*1,20 

420 WRITE! 6,403 ) TCENT ( I) , ( TCOUT ( I , J I , J* 1 , 13 ) 

401 F0RMAT!IH1,46H THERMAL CONDUCTIVITY * INTERNATIONAL BOOK 

402 FORMAT! 1H0,6HP- BARS, 1 3 ( 2 X, F7 .0 ) / 1H0, 3HT-F ) 

403 FORMAT ( 1H0 ,F5#0, 13F9. 2 I 
WRITE! 6,500 ) 

KU*3 

CO 450 J* 1 , 1 3 


PIN * PSVS!J) 

00 450 1*1,31 

TIN * TFVSm+460.0 

KR*0 

CALL WASP (1,16, TIN, PIN, D,M,KR) 

IF(KR.EQ.l) GO TO 450 

TC0UT2 ! I , J ) * K * 1. E 3 * 3600. 

450 CONTINUE 
WRITE!6,451I 

WRITE! 6,352) (PSVS! 11,1*1,13) 

CO 455 1*1,31 

455 WRITE! 6,353) TFVS ( I) , ( TCOUT 2 ( I , J ) , J* 1, 13) 

451 FORMAT! 1H1,20X,47H THERMAL CONDUCT I VI TY-ASME TABLE NO. 11 PAGE 281) 


SURFACE TENSION AND LAPLACE CONSTANT INTERNATIONAL BOOK 
PAGE 172 - TABLE 7 


WR I TE ( 6 ,481 ) 
KU*1 
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238 

r/ 7 


n n o 


DU 480 J« 1 » 16 
TI =• TCENT( Jl+273.15 
CALL WA$P(1*32»TI»PI » D t H f l ) 
NRITE(6*482) TCENT< J l »S IGMA» ALC 
F0RPAT(1W1#20X#75H INTERNATIONAL BOOK 
1 ION AND LAPLACE CONSTANT /1H0 # 24H 

A02 F0RPATI1H * F7.0* F8.2, F8.3 > 


460 

481 


— TABLE NO. 7* 
T-C OYNE/CN 


SURFACE 

MH 


TENS 

/ 5 


END 


TEST PROGRAM 


500 FORPAT(lHl) 
STOP 
END 


239 
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241 

242 

243 
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249 

250 
25 ’. 
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COMPARISON POINTS FOR ASME TABLE 3, REF. 1 
(VOLUME, ENTHALPY, ENTROPY FOR PRESSURE LISTED) 
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COMPARISON POINTS FOR ASME TABLE 10, REF. 
(VISCOSITY, 10' 5 slugs/ ( f t ) ( sec ) ) 
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COMPARISON POINTS FOR TABLE 6a, REF. 
(THERMAL CONDUCTIVITY, 10' 2 W/(m)(K)) 
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COMPARISON POINTS FOR ASME TABLE 11, REF. 
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COMPARISON POINTS FOR TABLE 7, REF. 2 
(SURFACE TENSION AND LAPLACE CONSTANT) 
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APPENDIX G 


METASTABLE SUBROUTINE (PM ETAS) 

Although property measurements for other than stable states are very difficult to 
make, metastable states are of interest in heat-transfer and fluid-flow calculations. 
The fundamental equation 


= ^ 0 (T) + RT[ln p + pQ(p, t)] 

represents a continuum of single-phase states between the saturated liquid and saturated 
vapor states which can be classified as either metastable or unstable (as for the Van der 
Waals equation). Consequently, properties of the superheated liquid and supersaturated 
vapor can be determined. It is pointed out in reference 3 that between 300° C and the 
critical temperature, the nonstable states as determined by the fundamental equation 
have a single maximum and minimum. At lower temperatures, more than one pair of 
extremum exists for which the authors of reference 3 attach no significance. 

The subroutine PMETAS (KU, T, D, P, KR) is provided to illustrate to the user of 
WASP how the metastable and unstable states can be determined. Given a density D and 
a temperature T (KU and KR have their usual meanings) the pressure P is returned 

P = pRT jl + pQ(p, t) + P 2 8< ^ X ^j 

The user can then formulate a locus of maximum and minimum points, as for the Van 
der Waals equation, and determine if the point is stable, metastable, or unstable. (Note 
that PMETAS will also return stable points provided D and T represent a stable point. ) 
Examples of the metastable and unstable loci are given as figures 14. 
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$1 

c 

c 

c 

c 

c 

c 

c 

c 

c 


BFTC PMETA 

SUBROUTINE PMET A S ( KU, T , C , P , KR ) 

4/ 10 _ 72 


THIS ROUTINE CALCULATES PRESSURE FOR ANY T INPUT ANO D INPUT. 

THIS ROUTINE DOES NOT DEFINE A REGION AND IS NOT 
CALLED BY -HASP-. THE USER CALLS IT DIRECTLY AND 
IT CAN BE USED IN THE METASTABLE STATE. 


C 

C 

C 

C 


COMMCN /CONV3/PCCNV (5) 

COMMON /TPARAM/TS 

COMMON /CCNSTS/ T AUC, RHCA , RHOB, TAUA , E ,R 

COMMON /CHECKS/ CCH l « CCH2 » PC HI , PCH2, PCH3 , TCH1 t TCH2 t TCH3 * DST , TST . 
1HSCH1 ,HSCH2 
COMMON / 1 ERROR/ I ROUT 
DS= CCHECK ( KU*D ) 

TS= TCHECK ( KU f K R , T ) 

IRCLT= 2 


ONE EQUATION FOR ALL REGIONS 

CALL QMUST ( DS ) 

CALL QMUST2 ( TS) 

PS* 1C00.*R*DS/TS*( l.+DS*(QCALC(TS)+DS*QDT(DS,TSn ) 
90 P= PS+PCONV ( KU) 

RETURN 

END 
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APPENDIX H 


THERMODYNAMIC RELATIONS AND DERIVATIVES 

The symbols C , C y , H, P, R, S, T, and p have the same meaning as defined 
elsewhere in this report. The other symbols used exclusively in this appendix are de- 
fined as follows: 


A = E - TS 
E 

F = H - TS 


K 

V 


Helmholtz free energy or work content 
internal energy 

Gibbs free energy or free energy 
equilibrium constant 
specific volume 


To illustrate the facility of the partial derivatives, Roder and Weber (ref. 17) give 
five which are useful to engineers: 


Specific heat input: 



Energy derivative: 



Isothermal bulk modulus: 
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Volume expansivity: 



The background material necessary to derive these and other parameters as the 
Joule- Thom s< m coefficient 


£ - 1 

X J 

can be found in most thermodynamic texts. 

WASP provides the partial derivatives (3P/3p) T and (3P/3T) . With the aid of the 
following thermodynamic derivatives and the Bridgeman Tables, any thermodynamic 
parameter can be found. The following thermodynamic tables were taken from refer- 
ence 18. 

Differential energy formulas: 


PC T 


3P 

T xclT 

IP (dP 

dp 


dE = T dS - P dV 
dH = T dS + V dP 


dA = - S dT - P dV 


dF = - S dT + V dP 


Maxwell relations: 
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Energy-function derivatives: 




Heat-capacity relations: 



92 



Effect of P or V on H or E: 



Temperature effect on AF/T = - R In K: 




\ T / 

D 3 In K AH 

3T 

3T “'_2 

. 

P 


Partial molal quantities, where Y is any extensive quantity: 


y i = 


_3Y 

L0n, 


P, T, n2, 


Y = n lYl + n 2 y 2 + . 


X 




+ 


. =0 


frj\ 9 2 Y 

3n. / 3n. 3n. 

^ J / 1 ] 




The so-called Bridgeman Tables are summarized as follows: 
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TABLE I. - OPERATIONS SHEET FOR SUBROUTINE WASP 3 


COMMON/PROPTY/KU, DL, DV, HL, HV, S, SL, SV, CV, CVL, CVV, CP, CPL, CPV, GAMMA, GAMMAL, GAMMAV, 
CL, CVP, MU, MUL, MUV, K, KL, KV, SIGMA, EXCL, EXCV, EXCESK 
REAL MU, MUL, MUV, K, KL, KV 


CALL WASP (KS, KP, T, P, D, H, KR) 


1 I 
I I 
I I 
I I 
I I 


\\\\ 

\ \ N \ 

\ N \ X 
\ N \ \ 
\ \ ' 


\ 


\ \ 


\ \ 
\' 


\ 


\ 


V 



Liquid 


1 


Region 

KR-0 Unknown, check KR returned 
KR^l Saturation 
KR = 2 Liquid 
KR^3 Gas and/or fluid 
Enthalpy, J/g 

o 

Density, g/cm 

9 

Pressure, MN/m 

* Temperature, K 

Thermodynamic and transport properties 15 

Only P, p, and T returned 
Enthalpy, J/g; (H), (HL), (HV) 

Entropy, J/(g)(K); (S), (SL), (SV) 

Specific heat at constant volume, J/(g)(K); (CV), (CVL), (CVV) 
Specific heat at constant pressure, J/(g)(K); (CP), (CPL), (CPV) 
Ratio of specific heats, C p /C y ; (GAMMA), (GAMMAL), (GAMMAV) 
Sonic velocity, cm/sec; (C), (CL), (CVP) 

Dynamic viscosity, g/(cm)(sec); (MU), (MUL), (MUV) 

Thermal conductivity, W/(cm)(K); (K), (KL), (KV) 

Surface tension, dyne/cm; (SIGMA) 

Input specification of independent properties 


IVapor 


KP-0 


KP= 1 

H 

KP = 2 

S 

KP=4 

c, 


C , 


Y 


c 

KP=8 

M 

KP-16 

k 

KP-32 

u 


Entropy, S 


KS=1 

p - f(T, P); given T, P find p 

KS-2 

P = f(T,p); given T, p find P 

KS=3 

T = f(P ,p); given P,p find T 

KS = 4 

T, p = f(P, H); given P, H find 

KS = 5 

T, p = f(P,S); given P, S find 


Notes: 


1. The units indicator, KU, must be set such that 1 < KU < 5 or no valid property values can be determined. See table II. 

2. Reset KR * 1 for each call to WASP to be assured of nonsaturation calculations (unless T = T and P = P ). 

3. Sample problem: sat sat ’ 

COMMON/PROPTY/KU, etc. (as above) 

REAL MU, etc. {as above) 

KU-1 


KR-0 
T-773. 0 
D=0. 178 


Call WASP (2, 31, T, P, D, H, KR) 

WASP will return P = 40. MN/m^, KR = 3, H - 2902. 4, and the following values in COMMON: S = 5. 4689, CV = 2 4503 
CP = 5. 7893, GAMMA = 2. 363, C = 57415. 4, MU = 0. 3682x10 K = 0. 1534xl0~ 2 . 

KP input is ^ KP options if more than one property is requested. For example, if enthalpy and entropy are desired, set 
KP equal to 3. 
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TABLE II. - UNITS SPECIFICATION 


Physical quantity 

Units specification 

KU=1 

u 

io 

KU=3 

Temperature 

Density 

Pressure 

Enthalpy 

Entropy, specific heat 
Sonic velocity 
Dynamic viscosity 
Thermal conductivity 
Surface tension 

K 

g/cm 3 

MN/m 2 

joule/g 

joule/(g)(K) 

cm/sec 

g/(cm)(sec) 

joule/(cm)(sec)(K) 

dyne/cm 

K 

g/cm 3 

atmospheres 

joule/g 

joule/(g)(K) 

cm/sec 

g/(cm)(sec) 

joule/(cm)(sec)(K) 

dyne/cm 

°R 

lbm/ft 3 
psia 
Btu/lbm 
Btu/(lbm)(°R) 
ft/ sec 

lbm/(ft)(sec) 

Btu/(ft)(sec)(°R) 

Ibf/ft 


a KU=4, 5 permit the user to work in other units; however, the proper conversions 
must be entered into BLOCK DATA. To add special set of units for KU=4 or 


KU=5: 

(1) User's program must contain following COMMON 
/CONVl/DCONV(5) 


or modify BLOCK DATA subprogram directly 
(See appendix E for listing of BLOCK DATA. ) 


/CONV2/TCONV(5) 

/CONV3/PCONV(5) 

/CONV4/SCONV(5) 

/CONV5/CCONV(5) 

/CONV6/HCONV(5) 

/CONV7/MCONV(5) 

/CONV8/KCONV(5) 

/CONV 9/STCONV(5) 

REAL MCONV, KCONVj 
(2) Store conversion factors in fourth and/or fifth position of each array such 
that (D in input unit desired)/DCONV(4) = g/cm 3 , etc. All conversion 
factors must change input to units of KU=1. For output then, 

(D in g/cm 3 )*DCONV(4) = (D in desired units). 




TABLE III. - PROGRAM ASSEMBLY 



User's data cards 
Other user subroutines 
Calls to subroutine WASP 

f KU=, KR= 

[REAL MU, etc. 

COMMON/PROPTY/, etc. 

Beginning of user's property subroutine 

User's main program 

Job control cards and WASP subroutines 


The subroutines in WASP may be loaded in any order with respect to the 
user’s program. To run successfully, there must appear in at least one 
user subroutine the following: 

(1) COMMON/PROPTY/, etc. 

REAL MU, etc. 

(2) KU=1 (or 2, 3, 4, 5) 

(3) Other input variable specifications 

(4) Calls to subroutine WASP 

The COMMON/PROPTY/, of course, must be in the main program or 
subroutine where the user expects answers to be returned from WASP. 

It could be in several or all user subroutines. 
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TABLE IV. - COEFFICIENTS OF Q- FUNCTION, i^-FUNCTION, AND VAPOR PRESSURE CURVE 


(a) Coefficients of Q-function 


i 

j 


1 

2 

3 

4 

5 

6 

7 


Coefficients, A^ 

1 

29. 492937 

-5. 1985860 

6.8335354 

-C 

>. 1564104 

-6. 3972405 

-3.9661401 

-0.69048554 

2 

-132. 13917 

7.7779182 

-26. 149751 

-0. 72546108 

26. 409282 

1.5.453061 

2.7407416 

3 

274. 64632 

-33. 301902 

65. 326396 

-9.2734289 

-47. 740374 

-29. 142470 

-5. 1028070 

4 

-360.93828 

-16. 254622 

-26. 181978 

4.3125840 

56.323130 

29. 568796 

3.9636085 

5 

342. 18431 

-177. 31074 

( 

) 

C 



) 

C 


( 

) 

6 

-244, 50042 

127. 48742 











7 

155. 18535 

137. 46153 











8 

5. 9728487 

155. 97836 









- 


9 

-410. 30848 

337. 31180 

-137.46618 

6.7874983 

136. 87317 

79.847970 

13.041253 

10 

-416. 05860 

-209. 88866 

-733. 96848 

10. 401717 

645.81880 

399. 17570 

71. 531353 


(b) Coefficients of 


i^Q-function 


i 

Coefficients C- 

1 

1855.3865 

2 

3. 278642 

3 

-.00037903 

4 

46. 174 

5 

-1. 02117 


(c) Coefficients of vapor 


pressure curve 


i 

Coefficients 

1 

2. 9304370 

2 

-2309. 5789 

3 

. 34522497X10' 1 

4 

13621289X10' 3 

5 

. 25878044X10" 6 

6 

24709 162xl0' 9 

7 

. 95937646X10' 13 
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TABLE V. - NECESSARY AND OPTIONAL ROUTINES 


(a) Necessary routines 


NAME (* indicates multiple entry) 

Reason 

BLOCK DATA 

Stores coefficients for the fundamental equation 

•(CHECK, TCHECK, PCHECK, DC HECK) 

Performs region and limit checks for all subroutines; con- 
verts user’s units to internal program units 

ROOT 

ROOTX 

SOLVE 

Mathematical routines used in all iterative solutions neces- 
sary to calculation of properties 

QCALC 

•(QMUST, QMUST2) 
QTD 

*(QDTA, QDT) 
Q2T2D 

•(Q2DTA, Q2DT) 

* (Q2D2TA, Q2D2T) 

Q-function and derivatives used in equation-of-state calcu- 
lations (See equations used by all KS and KP options. ) 

DENS 

PSSS 

*(DSF, DDSF) 

Used for KS=1 request and to determine region number for 
most other KS and KP options 


(b) Optional routines 


NAME (* indicates 
multiple entry) 

KS or KP option in- 
volved 

Statement numbers 
in subroutine WASP 

Additional conditions for removal 

PRESS 

KS=2 

20 

None 

TEMP 

TSS 

*(TSSF, DTSSF) 
*(TSF, DTSF) 

KS=3 (also KS = 4 and 
KS=5) 

30 

40 

45 

Must also remove TEMPPH, TEMPPS, 
TSHF, and TPSF 

TEMPPH 
•(TSHF, TPSF) 

KS=4 

40 

None 

TEMPPS 
*(TSHF, TPSF) 

KS=5 

45 

None 

ENTH 

KP=1 (also KS=4) 

60 

40 

Must also remove TEMPPH and TSHF 

ENT 

KP=2 

KS=5 

100 

45 

Must also remove TEMPPS and TPSF 

CPPRL 

KP=4 

130 to 140 

None 

vise 

KP=8 

160 to 170 

None 

THERM 

KP=16 

180 to 190 

None 

SURF | 

KP=32 

240 

None 
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Excess viscosity, 17 - i? lf g/(cm)($ec) 



310° C. 

9. 87 MN/nr 


Curve from c 



Excess thermal conductivity, A - A n , w/(cmi(KI 



> Liquid} Saturation locus 
Curve from program WASP 



Reduced density, p/p c 


Figure 1 . - Excess thermal conductivity as function of reduced density for region 603 K <T< 673 K 
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(b) 6P/dT>p as function of temperature. 

Figure 4. - Pressure and the derivative 6P/6TL as function of tempera- 
ture for selected isochores. 
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(b) Cp/Cy as function of temperature 

Figure 9. - Specific heat at constant volume and specific-heat ratio Cp/C v 
as function of temperature for selected isobars 
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(a) Thermal conductivity as function of temperature for selected isobars. 
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(bi Approximation to anomalous behavior of thermal conductivity for water 
in the near-critical thermodynamic state. 

Figure 12. - Thermal conductivity as function of temperature for selected 
isobars and total thermal conductivity as function of density ratio for 
selected isotherms. 
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Figure 13. - Surface tension and Laplace constant as function of temper- 
ature. 
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igure 14. - Pressure-volume relations for metastable region along selec- 
ted isotherms. 
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Figure 17. - Temperature differences as function of temperature. Compar- 
ison of WASP calculated PVT to International Skeleton Tables. 
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